diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a808f6bf..73d86219 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -132,6 +132,10 @@ jobs: args: -p spatialrust-features --features feature-iss - name: spatialrust-features-iss-meta args: -p spatialrust --features feature-iss + - name: spatialrust-features-iss-gpu + args: -p spatialrust-features --features feature-iss-gpu + - name: spatialrust-features-iss-gpu-meta + args: -p spatialrust --features feature-iss-gpu - name: spatialrust-features-normal-orient args: -p spatialrust-features --features feature-normal-orient - name: spatialrust-features-normal-orient-meta diff --git a/CHANGELOG.md b/CHANGELOG.md index c5ce8d39..b20de50d 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 `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; + incomplete/malformed graphs invalidate all output. CPU defaults, ordered + host compaction and attribute gathering are retained. GPU selection can + differ at numerical thresholds/ties; diagnostic eigen readback is explicit. + - 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 diff --git a/README.md b/README.md index ca4f7389..d64a36b0 100644 --- a/README.md +++ b/README.md @@ -661,9 +661,17 @@ 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; ISS/MLS remain CPU. +RANSAC and ordinary descriptor defaults are retained; MLS remains CPU. See [FPFH measurements](notes/2026-10-12_gpu_fpfh_descriptors.md). +ISS keypoints have an explicit `feature-iss-gpu` path: +`detector.detect_with_gpu(&runtime, &cloud, salient_capacity, non_max_capacity)` +runs both radius graphs, covariance/eigenvalue saliency and local maxima on GPU. +CPU tree construction, ordered flag compaction and schema/attribute gathering +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). + ### Library Load or save by file extension: diff --git a/crates/spatialrust-features/Cargo.toml b/crates/spatialrust-features/Cargo.toml index 5e215758..b488a36f 100644 --- a/crates/spatialrust-features/Cargo.toml +++ b/crates/spatialrust-features/Cargo.toml @@ -12,6 +12,7 @@ description = "Feature estimation for SpatialRust point clouds" default = ["feature-normal"] feature-normal = [] feature-iss = [] +feature-iss-gpu = ["feature-iss", "dep:spatialrust-gpu"] feature-normal-orient = [] feature-boundary = [] feature-normal-gpu = ["feature-normal", "dep:spatialrust-gpu"] diff --git a/crates/spatialrust-features/src/iss.rs b/crates/spatialrust-features/src/iss.rs index 0d9a7dd7..d03b969d 100644 --- a/crates/spatialrust-features/src/iss.rs +++ b/crates/spatialrust-features/src/iss.rs @@ -56,6 +56,16 @@ pub struct IssKeypointResult { pub keypoints: PointCloud, } +/// Explicit GPU ISS detection with complete input/flag/status payload receipts. +#[cfg(feature = "feature-iss-gpu")] +#[derive(Clone, Debug, PartialEq)] +pub struct GpuIssKeypointResult { + /// Ordered original indices and gathered keypoint cloud. + pub result: IssKeypointResult, + /// Tree/query upload and selection/status download; uniform controls excluded. + pub transfers: spatialrust_core::TransferStats, +} + /// Intrinsic Shape Signatures keypoint detector. #[derive(Clone, Copy, Debug, PartialEq)] pub struct IssKeypointDetector { @@ -84,6 +94,64 @@ impl IssKeypointDetector { self.detect_with_workers(input, parallel_worker_count(input.len())) } + /// Explicitly uploads input and computes radius graphs, covariance, + /// eigenvalue saliency and non-maximum suppression on the supplied GPU. + /// CPU tree construction, ordered flag compaction and attribute gathering + /// remain host operations; no automatic CPU fallback occurs. + /// Both capacities include self and must fit ALL radius matches. Overflow + /// rejects instead of silently using truncated rows. Positions/extents must + /// be finite, radii must have positive finite nonzero f32 squares, and gamma + /// values must be finite in [0,1]. These checks are stricter than CPU ISS. + /// GPU f32 covariance/eigen arithmetic can change selection near threshold + /// and NMS ties compared with the CPU's f64 arithmetic. + /// A nonzero evaluated local extent whose f32 square underflows the normal + /// range returns a numerical error; no approximate zero-variance fallback. + #[cfg(feature = "feature-iss-gpu")] + pub fn detect_with_gpu( + &self, + runtime: &spatialrust_gpu::WgpuRuntime, + input: &PointCloud, + salient_capacity: u32, + non_max_capacity: u32, + ) -> SpatialResult { + let radii = [self.config.salient_radius, self.config.non_max_radius]; + let gamma = [self.config.gamma_21, self.config.gamma_32]; + if radii.iter().any(|&r| !r.is_finite() || r <= 0.0 || !(r * r).is_finite() || r * r == 0.0) + || gamma.iter().any(|g| !g.is_finite() || !(0.0..=1.0).contains(g)) + { + return Err(SpatialError::InvalidArgument("GPU ISS requires positive finite nonzero radius squares and finite gamma values in [0,1]".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 ISS requires finite XYZ".to_owned())); + } + let tree = KdTree::from_slices(x, y, z); + let device_input = spatialrust_gpu::GpuIssInput::upload(runtime, &tree)?; + let neighborhoods = device_input.neighborhoods_resident( + runtime, + self.config.salient_radius, + self.config.non_max_radius, + salient_capacity, + non_max_capacity, + )?; + let output = device_input.keypoints_resident( + runtime, + &neighborhoods, + self.config.gamma_21, + self.config.gamma_32, + u32::try_from(self.config.min_neighbors).unwrap_or(u32::MAX), + )?; + let (indices, downloaded) = output.download(runtime)?; + let indices: Vec<_> = indices.into_iter().map(|i| i as usize).collect(); + let keypoints = gather_indices(input, &indices)?; + let mut transfers = device_input.upload_transfers(); + transfers.record( + spatialrust_core::TransferDirection::DeviceToHost, + downloaded.device_to_host_bytes(), + ); + Ok(GpuIssKeypointResult { result: IssKeypointResult { indices, keypoints }, transfers }) + } + fn detect_with_workers( &self, input: &PointCloud, diff --git a/crates/spatialrust-features/src/lib.rs b/crates/spatialrust-features/src/lib.rs index f1d18706..7eba5cde 100644 --- a/crates/spatialrust-features/src/lib.rs +++ b/crates/spatialrust-features/src/lib.rs @@ -30,6 +30,9 @@ pub use boundary::{BoundaryConfig, BoundaryDetector, BoundaryResult}; #[cfg(feature = "feature-iss")] pub use iss::{IssKeypointConfig, IssKeypointDetector, IssKeypointResult}; +#[cfg(feature = "feature-iss-gpu")] +pub use iss::GpuIssKeypointResult; + #[cfg(feature = "feature-normal-orient")] pub use orient::{orient_normals_consistent, NormalOrientationConfig}; diff --git a/crates/spatialrust-features/tests/fixtures/iss_gpu_reference.rs b/crates/spatialrust-features/tests/fixtures/iss_gpu_reference.rs new file mode 100644 index 00000000..509754c4 --- /dev/null +++ b/crates/spatialrust-features/tests/fixtures/iss_gpu_reference.rs @@ -0,0 +1,158 @@ +use spatialrust_core::{HasPositions3, PointCloud}; +use spatialrust_features::{IssKeypointConfig, IssKeypointResult}; +use spatialrust_gpu::{GpuIssInput, GpuIssNeighborhoods, WgpuRuntime}; +use spatialrust_math::{symmetric_eigen3, Mat3}; +use spatialrust_search::{KdTree, RadiusSearchIndex}; + +pub struct Audit { + pub indices: Vec, + pub differences: usize, + pub max_absolute_eigen_error: f64, +} + +fn cpu_scores(cloud: &PointCloud, config: IssKeypointConfig, tree: &KdTree) -> Vec<[f64; 4]> { + let (x, y, z) = cloud.positions3().unwrap(); + (0..cloud.len()) + .map(|i| { + let neighbors = tree.radius_search(x[i], y[i], z[i], config.salient_radius); + if neighbors.len() < 3 || neighbors.len() < config.min_neighbors { + return [0.0; 4]; + } + let count = neighbors.len() as f64; + let mut mean = [0.0; 3]; + for n in &neighbors { + for (m, value) in mean.iter_mut().zip([x[n.index], y[n.index], z[n.index]]) { + *m += f64::from(value); + } + } + for m in &mut mean { + *m /= count; + } + let mut c = [0.0; 6]; + for n in &neighbors { + let d = [ + f64::from(x[n.index]) - mean[0], + f64::from(y[n.index]) - mean[1], + f64::from(z[n.index]) - mean[2], + ]; + for (sum, value) in c.iter_mut().zip([ + d[0] * d[0], + d[1] * d[1], + d[2] * d[2], + d[0] * d[1], + d[0] * d[2], + d[1] * d[2], + ]) { + *sum += value; + } + } + for value in &mut c { + *value *= 1.0 / count; + } + let e = symmetric_eigen3(Mat3::from_rows( + [c[0], c[3], c[4]], + [c[3], c[1], c[5]], + [c[4], c[5], c[2]], + )) + .eigenvalues; + let eligible = e[2] > 0.0 + && e[1] > 0.0 + && e[1] / e[2] < f64::from(config.gamma_21) + && e[0] / e[1] < f64::from(config.gamma_32); + [e[0], e[1], e[2], if eligible { 1.0 } else { 0.0 }] + }) + .collect() +} + +// Bound eigenvalue errors against f64 math, independently reproduce device NMS +// from its downloaded scores and a CPU radius graph, and explain every CPU/GPU +// selection difference by a threshold or competing-saliency uncertainty band. +pub fn audit( + runtime: &WgpuRuntime, + cloud: &PointCloud, + config: IssKeypointConfig, + input: &GpuIssInput, + rows: &GpuIssNeighborhoods, + cpu: &IssKeypointResult, +) -> Audit { + let output = input + .keypoints_resident( + runtime, + rows, + config.gamma_21, + config.gamma_32, + u32::try_from(config.min_neighbors).unwrap_or(u32::MAX), + ) + .unwrap(); + let (scores, receipt) = output.download_scores(runtime).unwrap(); + assert_eq!( + receipt.device_to_host_bytes(), + if cloud.is_empty() { 0 } else { (cloud.len() * 16 + 32) as u64 } + ); + let indices = output.download(runtime).unwrap().0; + let tree = KdTree::from_point_cloud(cloud).unwrap(); + let reference = cpu_scores(cloud, config, &tree); + let tolerances: Vec<_> = + reference.iter().map(|row| 1e-12 + 5e-6 * (row[0] + row[1] + row[2]).abs()).collect(); + let mut max_absolute_eigen_error = 0.0_f64; + let mut ambiguous = vec![false; cloud.len()]; + for (i, (a, b)) in scores.iter().zip(&reference).enumerate() { + if b.iter().all(|&v| v == 0.0) { + assert_eq!(a[3], 0.0, "unsupported/zero covariance must be ineligible at row {i}"); + } + for bin in 0..3 { + let error = (f64::from(a[bin]) - b[bin]).abs(); + max_absolute_eigen_error = max_absolute_eigen_error.max(error); + assert!( + error <= tolerances[i], + "eigen row={i} bin={bin} gpu={} cpu={} error={error} bound={}", + a[bin], + b[bin], + tolerances[i] + ); + } + let boundary = b[2].abs() <= tolerances[i] + || b[1].abs() <= tolerances[i] + || (b[1] - f64::from(config.gamma_21) * b[2]).abs() + <= (1.0 + f64::from(config.gamma_21)) * tolerances[i] + || (b[0] - f64::from(config.gamma_32) * b[1]).abs() + <= (1.0 + f64::from(config.gamma_32)) * tolerances[i]; + ambiguous[i] = boundary; + assert!( + f64::from(a[3]) == b[3] || boundary, + "unexplained eligibility difference at row {i}" + ); + } + let (x, y, z) = cloud.positions3().unwrap(); + let mut from_scores = vec![]; + let mut differences = 0; + for i in 0..cloud.len() { + let neighbors = tree.radius_search(x[i], y[i], z[i], config.non_max_radius); + let selected = scores[i][3] != 0.0 + && neighbors.iter().all(|n| { + n.index == i || scores[n.index][3] == 0.0 || scores[i][0] >= scores[n.index][0] + }); + if selected { + from_scores.push(i as u32); + } + if selected != cpu.indices.binary_search(&i).is_ok() { + differences += 1; + assert!( + ambiguous[i] + || neighbors.iter().any(|n| n.index != i + && (reference[n.index][3] != 0.0 || scores[n.index][3] != 0.0) + && (ambiguous[n.index] + || (f64::from(reference[i][0] as f32) + - f64::from(reference[n.index][0] as f32)) + .abs() + <= tolerances[i] + tolerances[n.index])), + "unexplained NMS difference at row {i}" + ); + } + } + assert_eq!( + indices, from_scores, + "device NMS must match independent selection from its own scores" + ); + Audit { indices, differences, max_absolute_eigen_error } +} diff --git a/crates/spatialrust-features/tests/gpu_iss.rs b/crates/spatialrust-features/tests/gpu_iss.rs new file mode 100644 index 00000000..e272dd31 --- /dev/null +++ b/crates/spatialrust-features/tests/gpu_iss.rs @@ -0,0 +1,289 @@ +#![cfg(all(feature = "feature-iss-gpu", not(target_arch = "wasm32")))] + +use spatialrust_core::{ + DType, FieldSemantic, PointBuffer, PointBufferSet, PointCloud, PointField, PointSchema, + SpatialMetadata, StandardSchemas, Timestamp, +}; +use spatialrust_features::{IssKeypointConfig, IssKeypointDetector, IssKeypointResult}; +use spatialrust_gpu::{GpuIssInput, WgpuRuntime}; + +#[path = "fixtures/iss_gpu_reference.rs"] +mod math_reference; +pub use math_reference::audit; + +pub fn fixture(n: usize, kind: 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 % 128) as f32 * 0.05 - 3.2; + let y = (j / 128) as f32 * 0.05 - 1.6; + let z = if kind == 0 { + 0.0 + } else { + 0.015 * (x * x + y * y) + ((j * 37 % 101) as f32 / 100.0 - 0.5) * 0.08 + }; + 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: &IssKeypointResult, expected: &IssKeypointResult) { + assert_eq!(actual.indices, expected.indices); + assert_eq!(actual.keypoints, expected.keypoints); +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_iss_compares_selection_and_receipts_with_current_cpu() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let mut decisions = 0; + let mut differences = 0; + let mut max_error = 0.0_f64; + for n in [0, 1, 2, 17, 63, 64, 65, 257, 4097] { + for kind in 0..3 { + let cloud = fixture(n, kind); + for radius in [0.075, 0.16] { + let detector = + IssKeypointDetector::new(IssKeypointConfig::with_radii(radius, radius * 0.75)); + let expected = detector.detect(&cloud).unwrap(); + let actual = detector.detect_with_gpu(&runtime, &cloud, 128, 128).unwrap(); + let tree = spatialrust_search::KdTree::from_point_cloud(&cloud).unwrap(); + let input = GpuIssInput::upload(&runtime, &tree).unwrap(); + let rows = input + .neighborhoods_resident(&runtime, radius, radius * 0.75, 128, 128) + .unwrap(); + let checked = audit(&runtime, &cloud, detector.config(), &input, &rows, &expected); + assert_eq!( + actual.result.indices, + checked.indices.iter().map(|&i| i as usize).collect::>() + ); + differences += checked.differences; + max_error = max_error.max(checked.max_absolute_eigen_error); + eprintln!("GPU_ISS_ACCURACY n={n} kind={kind} radius={radius} cpu={} gpu={} differences={} max_abs_eigen_error={}",expected.indices.len(),actual.result.indices.len(),checked.differences,checked.max_absolute_eigen_error); + decisions += n; + let nodes = spatialrust_search::KdTree::from_point_cloud(&cloud) + .unwrap() + .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 * 4 + 32) as u64 } + ); + } + } + } + eprintln!( + "GPU_ISS audited {decisions} decisions and {} eigenvalues; {differences} CPU/GPU selection differences confined to bounded numerical boundaries; max_abs_eigen_error={max_error}",decisions*3 + ); +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_iss_validates_configuration_capacity_degeneracy_and_huge_radii() { + 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] { + for config in [ + IssKeypointConfig::with_radii(radius, 0.12), + IssKeypointConfig::with_radii(0.16, radius), + ] { + assert!(IssKeypointDetector::new(config) + .detect_with_gpu(&runtime, &cloud, 128, 128) + .is_err()); + } + } + for gamma in [-0.1, 1.1, f32::NAN, f32::INFINITY] { + for config in [ + IssKeypointConfig { gamma_21: gamma, ..Default::default() }, + IssKeypointConfig { gamma_32: gamma, ..Default::default() }, + ] { + assert!(IssKeypointDetector::new(config) + .detect_with_gpu(&runtime, &cloud, 128, 128) + .is_err()); + } + } + let detector = IssKeypointDetector::new(IssKeypointConfig::with_radii(0.16, 0.12)); + for (salient, non_max) in [(0, 128), (1, 128), (128, 0), (128, 1)] { + assert!(detector.detect_with_gpu(&runtime, &cloud, salient, non_max).is_err()); + } + for columns in [ + [vec![0.0; 33], vec![0.0; 33], vec![0.0; 33]], + [(0..33).map(|i| i as f32 * 0.01).collect(), vec![0.0; 33], vec![0.0; 33]], + [vec![f32::NAN], vec![0.0], vec![0.0]], + [vec![f32::MAX, -f32::MAX], vec![0.0; 2], vec![0.0; 2]], + ] { + let cloud = from_columns(columns); + let result = detector.detect_with_gpu(&runtime, &cloud, 128, 128); + if cloud.len() == 33 { + assert!(result.unwrap().result.indices.is_empty()); + } else { + assert!(result.is_err()); + } + } + let sparse = fixture(17, 1); + for config in [ + IssKeypointConfig { min_neighbors: usize::MAX, ..Default::default() }, + IssKeypointConfig { gamma_21: 0.0, ..Default::default() }, + ] { + assert!(IssKeypointDetector::new(config) + .detect_with_gpu(&runtime, &sparse, 128, 128) + .unwrap() + .result + .indices + .is_empty()); + } + let huge = IssKeypointDetector::new(IssKeypointConfig::with_radii(1e19, 1e19)); + compare( + &huge.detect_with_gpu(&runtime, &sparse, 32, 32).unwrap().result, + &huge.detect(&sparse).unwrap(), + ); + // Local variance too small for a normal f32 square is rejected explicitly. + let tiny = from_columns([vec![0.0, 1e-22, 2e-22], vec![0.0; 3], vec![0.0; 3]]); + assert!(IssKeypointDetector::new(IssKeypointConfig { + min_neighbors: 3, + ..IssKeypointConfig::with_radii(1.0, 1.0) + }) + .detect_with_gpu(&runtime, &tiny, 3, 3) + .is_err()); +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_iss_preserves_all_scalar_attribute_bits_semantics_and_metadata() { + let source = fixture(4097, 1); + let mut schema = PointSchema::new(); + let mut buffers = PointBufferSet::new(); + for (field, name) in source.schema().fields().iter().zip(["px", "py", "pz"]) { + schema = schema.with_field(PointField::scalar(name, field.semantic, field.dtype)); + buffers.insert(name, source.field(&field.name).unwrap().clone()); + } + for (name, dtype, buffer) in [ + ( + "float", + DType::F32, + PointBuffer::F32( + (0..source.len()) + .map(|i| f32::from_bits(if i % 2 == 0 { 0x7fc00037 } else { 0x80000000 })) + .collect(), + ), + ), + ( + "double", + DType::F64, + PointBuffer::F64( + (0..source.len()) + .map(|i| { + f64::from_bits(if i % 2 == 0 { + 0x7ff8000000000037 + } else { + 0x8000000000000000 + }) + }) + .collect(), + ), + ), + ("byte", DType::U8, PointBuffer::U8((0..source.len()).map(|i| i as u8).collect())), + ("short", DType::U16, PointBuffer::U16((0..source.len()).map(|i| i as u16).collect())), + ("index", DType::U32, PointBuffer::U32((0..source.len()).map(|i| i as u32).collect())), + ("signed", DType::I32, PointBuffer::I32((0..source.len()).map(|i| -(i as i32)).collect())), + ] { + schema = schema.with_field(PointField::scalar(name, FieldSemantic::Unknown, dtype)); + buffers.insert(name, buffer); + } + let input = PointCloud::try_from_parts( + schema, + buffers, + SpatialMetadata::new("sensor", Timestamp::from_nanos(37)), + ) + .unwrap(); + let runtime = WgpuRuntime::new_headless().unwrap(); + let detector = IssKeypointDetector::new(IssKeypointConfig::with_radii(0.16, 0.12)); + let actual = detector.detect_with_gpu(&runtime, &input, 128, 128).unwrap().result; + let expected = detector.detect(&input).unwrap(); + assert_eq!(actual.indices, expected.indices); + assert!(!actual.indices.is_empty()); + assert_eq!(actual.keypoints.schema(), input.schema()); + assert_eq!(actual.keypoints.metadata(), input.metadata()); + for field in input.schema().fields() { + match ( + actual.keypoints.field(&field.name).unwrap(), + expected.keypoints.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), + } + } +} + +#[test] +#[ignore = "requires a native wgpu adapter; run explicitly with --include-ignored"] +fn gpu_iss_reuses_owned_resident_outputs_and_checks_runtime_and_frame() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let other = WgpuRuntime::new_headless().unwrap(); + let cloud = fixture(257, 1); + let detector = IssKeypointDetector::new(IssKeypointConfig::with_radii(0.16, 0.12)); + let tree = spatialrust_search::KdTree::from_point_cloud(&cloud).unwrap(); + let input = GpuIssInput::upload(&runtime, &tree).unwrap(); + let different = GpuIssInput::upload(&runtime, &tree).unwrap(); + let rows = input.neighborhoods_resident(&runtime, 0.16, 0.12, 128, 128).unwrap(); + let expected = audit( + &runtime, + &cloud, + detector.config(), + &input, + &rows, + &detector.detect(&cloud).unwrap(), + ) + .indices; + assert!(input.neighborhoods_resident(&other, 0.16, 0.12, 128, 128).is_err()); + assert!(different.keypoints_resident(&runtime, &rows, 0.975, 0.975, 5).is_err()); + assert!(input.keypoints_resident(&other, &rows, 0.975, 0.975, 5).is_err()); + let retained = input.keypoints_resident(&runtime, &rows, 0.975, 0.975, 5).unwrap(); + std::thread::scope(|scope| { + let calls: Vec<_> = (0..4) + .map(|_| { + scope.spawn(|| { + for _ in 0..2 { + assert_eq!( + input + .keypoints_resident(&runtime, &rows, 0.975, 0.975, 5) + .unwrap() + .download(&runtime) + .unwrap() + .0, + expected + ); + } + }) + }) + .collect(); + for call in calls { + call.join().unwrap(); + } + }); + drop(rows); + drop(input); + drop(tree); + assert_eq!(retained.download(&runtime).unwrap().0, expected); + assert!(retained.download(&other).is_err()); +} diff --git a/crates/spatialrust-features/tests/gpu_iss_benchmark.rs b/crates/spatialrust-features/tests/gpu_iss_benchmark.rs new file mode 100644 index 00000000..2638d081 --- /dev/null +++ b/crates/spatialrust-features/tests/gpu_iss_benchmark.rs @@ -0,0 +1,106 @@ +#![cfg(all(feature = "feature-iss-gpu", not(target_arch = "wasm32")))] + +#[path = "gpu_iss.rs"] +mod reference; + +use spatialrust_features::{IssKeypointConfig, IssKeypointDetector}; +use spatialrust_gpu::{GpuIssInput, WgpuRuntime}; +use std::time::Instant; + +#[test] +#[ignore = "completed ISS benchmark; run this exact test alone in release mode"] +fn completed_iss_calls() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let info = runtime.adapter_info(); + println!("GPU_ISS_ENV backend={} device_type={}", info.backend, info.device_type); + for n in [257, 4097, 16384, 65537] { + for kind in 0..3 { + let radius = if kind == 0 { 0.11 } else { 0.16 }; + let salient_capacity = if kind == 2 { 96 } else { 64 }; + let non_max_capacity = 64; + let cloud = reference::fixture(n, kind); + let config = IssKeypointConfig::with_radii(radius, radius * 0.75); + let detector = IssKeypointDetector::new(config); + let cpu_expected = detector.detect(&cloud).unwrap(); + let tree = spatialrust_search::KdTree::from_point_cloud(&cloud).unwrap(); + let input = GpuIssInput::upload(&runtime, &tree).unwrap(); + let rows = input + .neighborhoods_resident( + &runtime, + radius, + radius * 0.75, + salient_capacity, + non_max_capacity, + ) + .unwrap(); + let checked = reference::audit(&runtime, &cloud, config, &input, &rows, &cpu_expected); + let gpu_expected = detector + .detect_with_gpu(&runtime, &cloud, salient_capacity, non_max_capacity) + .unwrap() + .result; + assert_eq!( + gpu_expected.indices, + checked.indices.iter().map(|&i| i as usize).collect::>() + ); + let cpu = || detector.detect(&cloud).unwrap(); + let gpu = || { + detector + .detect_with_gpu(&runtime, &cloud, salient_capacity, non_max_capacity) + .unwrap() + .result + }; + for _ in 0..3 { + reference::compare(&cpu(), &cpu_expected); + reference::compare(&gpu(), &gpu_expected); + assert_eq!( + input + .keypoints_resident( + &runtime, + &rows, + config.gamma_21, + config.gamma_32, + config.min_neighbors as u32 + ) + .unwrap() + .download(&runtime) + .unwrap() + .0, + checked.indices + ); + } + 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 + .keypoints_resident( + &runtime, + &rows, + config.gamma_21, + config.gamma_32, + config.min_neighbors as u32, + ) + .unwrap(); + runtime.wait_idle(); + resident_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output.download(&runtime).unwrap().0, checked.indices); + } else { + let result = if stage == 0 { cpu() } else { gpu() }; + let elapsed = start.elapsed().as_secs_f64() * 1000.0; + if stage == 0 { + cpu_ms.push(elapsed); + reference::compare(&result, &cpu_expected); + } else { + gpu_ms.push(elapsed); + reference::compare(&result, &gpu_expected); + } + } + } + } + let workers = spatialrust_search::parallel_worker_count(n); + println!("GPU_ISS_BENCH n={n} kind={kind} workers={workers} salient_capacity={salient_capacity} non_max_capacity={non_max_capacity} cpu_selected={} gpu_selected={} selection_differences={} max_abs_eigen_error={} cpu_ms={cpu_ms:?} gpu_ms={gpu_ms:?} resident_ms={resident_ms:?}",cpu_expected.indices.len(),gpu_expected.indices.len(),checked.differences,checked.max_absolute_eigen_error); + } + } +} diff --git a/crates/spatialrust-gpu/src/kernels/iss.rs b/crates/spatialrust-gpu/src/kernels/iss.rs new file mode 100644 index 00000000..f19fd6f1 --- /dev/null +++ b/crates/spatialrust-gpu/src/kernels/iss.rs @@ -0,0 +1,469 @@ +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::{GpuKdTree, GpuNeighborhoods, GpuQueryBatch, WgpuRuntime}; + +#[repr(C)] +#[derive(Clone, Copy, Pod, Zeroable)] +struct Params { + count: u32, + salient_capacity: u32, + non_max_capacity: u32, + min_neighbors: u32, + gamma_21: f32, + gamma_32: f32, + radius: f32, + non_max_radius_squared: f32, +} + +fn invalid() -> SpatialError { + SpatialError::InvalidArgument("GPU ISS requires matching input/runtime ownership, finite XYZ/extents, positive finite nonzero radius squares, finite gamma values in [0,1] 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(16).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 < 32 + || 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 ISS positions/tree/queries, reusable on one runtime. +/// CPU tree construction is external; subsequent search and selection run on GPU. +pub struct GpuIssInput { + tree: GpuKdTree, + queries: GpuQueryBatch, + count: u32, + identity: Arc<()>, + input_identity: Arc<()>, + transfers: TransferStats, +} + +impl GpuIssInput { + /// 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 + } + + /// Computes both reusable radius graphs, including self, without readback + /// or a host wait. Every match is counted; selection rejects overflow. + pub fn neighborhoods_resident( + &self, + runtime: &WgpuRuntime, + salient_radius: f32, + non_max_radius: f32, + salient_capacity: u32, + non_max_capacity: u32, + ) -> SpatialResult { + if !Arc::ptr_eq(&self.identity, &runtime.profiling_identity()) + || !valid_radius(salient_radius) + || !valid_radius(non_max_radius) + { + return Err(invalid()); + } + let salient = self.tree.radius_neighbors_resident( + runtime, + &self.queries, + salient_radius, + salient_capacity, + )?; + let non_max = self.tree.radius_neighbors_resident( + runtime, + &self.queries, + non_max_radius, + non_max_capacity, + )?; + Ok(GpuIssNeighborhoods { + salient, + non_max, + salient_radius, + non_max_radius, + input_identity: Arc::clone(&self.input_identity), + }) + } + + /// Runs covariance/eigenvalue saliency then local-maximum selection in two + /// passes in one submission. Both graphs must be complete, even for rejected + /// candidates. Gamma values must be finite in [0,1]; min_neighbors includes + /// self. Equal computed maxima all survive. No payload transfer or host wait. + /// GPU f32 arithmetic can change threshold/NMS decisions versus CPU f64. + /// Nonzero local extents whose f32 square is subnormal/zero invalidate the + /// result with numerical status rather than evaluating underflowed covariance. + /// Bad graph status invalidates and zeros ALL scores/flags; GPU consumers + /// must check status and explicit download rejects invalid output. + pub fn keypoints_resident( + &self, + runtime: &WgpuRuntime, + neighbors: &GpuIssNeighborhoods, + gamma_21: f32, + gamma_32: f32, + min_neighbors: u32, + ) -> SpatialResult { + if !Arc::ptr_eq(&self.identity, &runtime.profiling_identity()) + || !Arc::ptr_eq(&self.input_identity, &neighbors.input_identity) + || [gamma_21, gamma_32].iter().any(|g| !g.is_finite() || !(0.0..=1.0).contains(g)) + { + 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 = GpuIssKeypoints { + scores: create("iss-eigenvalues", (u64::from(self.count) * 16).max(16)), + flags: create("iss-selection", (u64::from(self.count) * 4).max(4)), + status: create("iss-status", 32), + count: self.count, + identity: Arc::clone(&self.identity), + }; + if self.count == 0 { + return Ok(result); + } + let cache = runtime.iss_pipelines.get_or_init(|| IssPipelines::new(runtime.device())); + let params = Params { + count: self.count, + salient_capacity: neighbors.salient.capacity(), + non_max_capacity: neighbors.non_max.capacity(), + min_neighbors, + gamma_21, + gamma_32, + radius: neighbors.salient_radius, + non_max_radius_squared: neighbors.non_max_radius * neighbors.non_max_radius, + }; + let uniform = runtime.device().create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some("iss-params"), + contents: bytemuck::bytes_of(¶ms), + usage: wgpu::BufferUsages::UNIFORM, + }); + let buffers = [ + &uniform, + self.queries.buffer().unwrap(), + neighbors.salient.buffer(), + neighbors.non_max.buffer(), + &result.scores, + &result.flags, + &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("iss-bindings"), + layout: &cache.layout, + entries: &entries, + }); + let mut encoder = runtime.device().create_command_encoder(&Default::default()); + for pipeline in [&cache.saliency, &cache.selection] { + let mut pass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("iss-compute"), + timestamp_writes: None, + }); + pass.set_pipeline(pipeline); + pass.set_bind_group(0, &group, &[]); + pass.dispatch_workgroups(self.count.div_ceil(64), 1, 1); + } + runtime.queue().submit(Some(encoder.finish())); + Ok(result) + } +} + +/// Two radius graphs tied to exactly one immutable uploaded ISS input. +pub struct GpuIssNeighborhoods { + salient: GpuNeighborhoods, + non_max: GpuNeighborhoods, + salient_radius: f32, + non_max_radius: f32, + input_identity: Arc<()>, +} +impl GpuIssNeighborhoods { + /// Covariance radius graph. Treat exposed storage as immutable. + #[must_use] + pub fn salient_rows(&self) -> &GpuNeighborhoods { + &self.salient + } + /// Non-maximum-suppression graph. Treat exposed storage as immutable. + #[must_use] + pub fn non_max_rows(&self) -> &GpuNeighborhoods { + &self.non_max + } +} + +/// Owned resident ISS selection flags, eigenvalue/eligibility rows and validity +/// status. No result storage is pooled. Custom consumers must inspect status. +pub struct GpuIssKeypoints { + scores: wgpu::Buffer, + flags: wgpu::Buffer, + status: wgpu::Buffer, + count: u32, + identity: Arc<()>, +} +impl GpuIssKeypoints { + /// Original-order u32 flags (0 or 1), with 4 unused bytes for empty input. + #[must_use] + pub fn buffer(&self) -> &wgpu::Buffer { + &self.flags + } + /// Original-order vec4: ascending eigenvalues and eligibility (0 or 1). + /// Rejected ratio tests still expose eigenvalues; insufficient support is + /// zero. Empty input has 16 unused bytes. Invalid graph output is all zero. + #[must_use] + pub fn scores_buffer(&self) -> &wgpu::Buffer { + &self.scores + } + /// Eight u32 words: bad salient rows, bad NMS rows, required salient/NMS + /// capacity, malformed rows, last bad row plus one, numerical errors, reserved. + /// All zero means valid output. Consumers must reject any nonzero status. + #[must_use] + pub fn status_buffer(&self) -> &wgpu::Buffer { + &self.status + } + /// Original point/flag count. + #[must_use] + pub const fn point_count(&self) -> u32 { + self.count + } + /// Explicit flag/status readback payload; scores are not downloaded. + #[must_use] + pub const fn readback_bytes(&self) -> u64 { + if self.count == 0 { + 0 + } else { + self.count as u64 * 4 + 32 + } + } + /// Explicitly downloads flags/status and returns ascending original indices. + /// Waits for completion; malformed/incomplete graphs reject before decoding. + /// Payload is 4 bytes per input point plus 32 status bytes, empty input zero. + /// CPU ordered flag compaction occurs here; no neighborhood/eigen readback. + pub fn download(&self, runtime: &WgpuRuntime) -> SpatialResult<(Vec, TransferStats)> { + self.readback(runtime, false, |words| { + if words.iter().any(|&v| v > 1) { + return Err(invalid()); + } + Ok(words + .iter() + .enumerate() + .filter_map(|(i, &flag)| (flag == 1).then_some(i as u32)) + .collect()) + }) + } + + /// Explicitly downloads eigenvalue/eligibility rows and status for diagnostics. + /// This separate opt-in transfer costs 16 bytes per point plus 32 status + /// bytes (empty input zero); the ordinary selection API does not call it. + /// Waits for completion and rejects invalid status/nonfinite or unsorted rows. + pub fn download_scores( + &self, + runtime: &WgpuRuntime, + ) -> SpatialResult<(Vec<[f32; 4]>, TransferStats)> { + self.readback(runtime, true, |words| { + let mut rows = Vec::with_capacity(self.count as usize); + for row in words.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()) + || values[0] > values[1] + || values[1] > values[2] + || (values[3] != 0.0 && values[3] != 1.0) + { + return Err(invalid()); + } + rows.push(values); + } + Ok(rows) + }) + } + + fn readback( + &self, + runtime: &WgpuRuntime, + scores: bool, + decode: impl FnOnce(&[u32]) -> SpatialResult, + ) -> SpatialResult<(T, TransferStats)> { + if !Arc::ptr_eq(&self.identity, &runtime.profiling_identity()) { + return Err(invalid()); + } + if self.count == 0 { + return Ok((T::default(), TransferStats::default())); + } + let bytes = u64::from(self.count) * if scores { 16 } else { 4 }; + let staging = runtime.device().create_buffer(&wgpu::BufferDescriptor { + label: Some("iss-readback"), + size: bytes + 32, + 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( + if scores { &self.scores } else { &self.flags }, + 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 payload_words = bytes as usize / 4; + let status = &words[payload_words..]; + let decoded = if status.iter().any(|&v| v != 0) { + Err(SpatialError::InvalidArgument(format!("GPU ISS rejected incomplete/malformed neighborhoods or numerical output; status={status:?}"))) + } else { + decode(&words[..payload_words]) + }; + drop(mapped); + staging.unmap(); + let mut transfers = TransferStats::default(); + transfers.record(TransferDirection::DeviceToHost, bytes + 32); + Ok((decoded?, transfers)) + } +} + +pub(crate) struct IssPipelines { + layout: wgpu::BindGroupLayout, + saliency: wgpu::ComputePipeline, + selection: wgpu::ComputePipeline, +} +impl IssPipelines { + 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 < 4 } + }, + has_dynamic_offset: false, + min_binding_size: None, + }, + count: None, + }) + .collect(); + let layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { + label: Some("iss-layout"), + entries: &entries, + }); + let pipeline_layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor { + label: Some("iss-layout"), + bind_group_layouts: &[&layout], + push_constant_ranges: &[], + }); + let shader = device.create_shader_module(wgpu::ShaderModuleDescriptor { + label: Some("iss-shader"), + source: wgpu::ShaderSource::Wgsl(include_str!("../shaders/iss.wgsl").into()), + }); + let create = |entry| { + device.create_compute_pipeline(&wgpu::ComputePipelineDescriptor { + label: Some("iss-pipeline"), + layout: Some(&pipeline_layout), + module: &shader, + entry_point: Some(entry), + compilation_options: Default::default(), + cache: None, + }) + }; + Self { saliency: create("build_saliency"), selection: create("select_keypoints"), layout } + } +} + +#[cfg(all(test, not(target_arch = "wasm32")))] +#[path = "iss_tests.rs"] +mod tests; diff --git a/crates/spatialrust-gpu/src/kernels/iss_tests.rs b/crates/spatialrust-gpu/src/kernels/iss_tests.rs new file mode 100644 index 00000000..eaf955e3 --- /dev/null +++ b/crates/spatialrust-gpu/src/kernels/iss_tests.rs @@ -0,0 +1,184 @@ +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 +} + +#[test] +fn iss_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) / 16) + .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 = 31, + _ => restricted.max_compute_workgroups_per_dimension = 0, + } + assert!(plan(1, &restricted).is_err()); + } +} + +#[test] +fn iss_validates_upload_configuration_and_empty_receipts() { + 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!(GpuIssInput::upload(&runtime, &KdTree::from_slices(&x, &zero, &zero)).is_err()); + } + let input = GpuIssInput::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, 1.0, 0, 0).is_err()); + assert!(input.neighborhoods_resident(&runtime, 1.0, radius, 0, 0).is_err()); + } + let rows = input.neighborhoods_resident(&runtime, 1.0, 1.0, 0, 0).unwrap(); + for gamma in [-0.1, 1.1, f32::NAN, f32::INFINITY] { + assert!(input.keypoints_resident(&runtime, &rows, gamma, 0.9, 0).is_err()); + assert!(input.keypoints_resident(&runtime, &rows, 0.9, gamma, 0).is_err()); + } + let output = input.keypoints_resident(&runtime, &rows, 0.9, 0.9, 0).unwrap(); + assert_eq!(output.point_count(), 0); + assert_eq!(output.readback_bytes(), 0); + assert_eq!(output.download(&runtime).unwrap(), (vec![], TransferStats::default())); + assert_eq!(output.download_scores(&runtime).unwrap(), (vec![], TransferStats::default())); +} + +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 incomplete_or_poisoned_graphs_zero_every_flag_and_score() { + 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 = GpuIssInput::upload(&runtime, &KdTree::from_slices(&x, &zero, &zero)).unwrap(); + for second in [false, true] { + let rows = if second { + input.neighborhoods_resident(&runtime, 0.2, 100.0, 3, 3).unwrap() + } else { + input.neighborhoods_resident(&runtime, 100.0, 0.01, 3, 1).unwrap() + }; + let output = input.keypoints_resident(&runtime, &rows, 0.975, 0.975, 3).unwrap(); + assert_eq!( + read_words(&runtime, output.status_buffer(), 32), + if second { [0, 65, 0, 65, 0, 65, 0, 0] } else { [65, 0, 65, 0, 0, 65, 0, 0] } + ); + assert!(read_words(&runtime, output.buffer(), 65 * 4).iter().all(|&v| v == 0)); + assert!(read_words(&runtime, output.scores_buffer(), 65 * 16).iter().all(|&v| v == 0)); + assert!(output.download(&runtime).is_err()); + assert!(output.download_scores(&runtime).is_err()); + let rows = input.neighborhoods_resident(&runtime, 0.2, 0.01, 3, 1).unwrap(); + poison( + &runtime, + if second { rows.non_max_rows().buffer() } else { rows.salient_rows().buffer() }, + ); + let output = input.keypoints_resident(&runtime, &rows, 0.975, 0.975, 3).unwrap(); + assert_eq!( + read_words(&runtime, output.status_buffer(), 32), + if second { [0, 1, 0, 1, 1, 1, 0, 0] } else { [1, 0, 2, 0, 1, 1, 0, 0] } + ); + assert!(read_words(&runtime, output.buffer(), 65 * 4).iter().all(|&v| v == 0)); + assert!(read_words(&runtime, output.scores_buffer(), 65 * 16).iter().all(|&v| v == 0)); + assert!(output.download(&runtime).is_err()); + } + let rows = input.neighborhoods_resident(&runtime, 0.2, 0.01, 3, 1).unwrap(); + let good = input.keypoints_resident(&runtime, &rows, 0.975, 0.975, 3).unwrap(); + assert_eq!(read_words(&runtime, good.status_buffer(), 32), [0; 8]); + assert!(good.download(&runtime).unwrap().0.is_empty()); +} + +#[test] +fn analytic_box_eigenvalues_and_equal_maxima_survive() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let (mut x, mut y, mut z) = (vec![], vec![], vec![]); + for px in [-0.06, 0.06] { + for py in [-0.04, 0.04] { + for pz in [-0.02, 0.02] { + x.push(px); + y.push(py); + z.push(pz); + } + } + } + x.push(100.0); + y.push(100.0); + z.push(100.0); + let tree = KdTree::from_slices(&x, &y, &z); + let input = GpuIssInput::upload(&runtime, &tree).unwrap(); + let rows = input.neighborhoods_resident(&runtime, 0.2, 1.0, 9, 9).unwrap(); + let output = input.keypoints_resident(&runtime, &rows, 0.975, 0.975, 5).unwrap(); + assert_eq!(output.download(&runtime).unwrap().0, (0..8).collect::>()); + let scores = read_words(&runtime, output.scores_buffer(), 9 * 16); + for row in scores[..32].chunks_exact(4) { + for (&bits, expected) in row[..3].iter().zip([0.0004, 0.0016, 0.0036]) { + assert!((f32::from_bits(bits) - expected).abs() < 2e-9); + } + assert_eq!(f32::from_bits(row[3]), 1.0); + } + assert_eq!(scores[32..], [0; 4]); + for (gamma, min_neighbors) in [(0.4, 5), (0.975, 9)] { + assert!(input + .keypoints_resident(&runtime, &rows, gamma, 0.975, min_neighbors) + .unwrap() + .download(&runtime) + .unwrap() + .0 + .is_empty()); + } +} diff --git a/crates/spatialrust-gpu/src/kernels/mod.rs b/crates/spatialrust-gpu/src/kernels/mod.rs index 5e8c3d22..5903f53a 100644 --- a/crates/spatialrust-gpu/src/kernels/mod.rs +++ b/crates/spatialrust-gpu/src/kernels/mod.rs @@ -51,6 +51,8 @@ mod voxel_gather; pub(crate) mod cluster_union; pub(crate) mod fpfh; +pub(crate) mod iss; +pub use iss::{GpuIssInput, GpuIssKeypoints, GpuIssNeighborhoods}; pub(crate) mod nearest; pub use fpfh::{GpuFpfhDescriptors, GpuFpfhInput, GpuFpfhNeighborhoods, GPU_FPFH_DESCRIPTOR_LEN}; pub use nearest::{ diff --git a/crates/spatialrust-gpu/src/kernels/neighborhoods.rs b/crates/spatialrust-gpu/src/kernels/neighborhoods.rs index 3226c3cc..61a35f65 100644 --- a/crates/spatialrust-gpu/src/kernels/neighborhoods.rs +++ b/crates/spatialrust-gpu/src/kernels/neighborhoods.rs @@ -26,6 +26,10 @@ pub struct GpuQueryBatch { } impl GpuQueryBatch { + pub(crate) fn buffer(&self) -> Option<&wgpu::Buffer> { + self.buffer.as_ref() + } + /// Uploads XYZ queries as padded 16-byte records. Finite XYZ and supported /// buffer/dispatch sizes are required; empty batches allocate no buffer. pub fn upload(runtime: &WgpuRuntime, queries: &[[f32; 3]]) -> SpatialResult { diff --git a/crates/spatialrust-gpu/src/lib.rs b/crates/spatialrust-gpu/src/lib.rs index 737ff219..5346610f 100644 --- a/crates/spatialrust-gpu/src/lib.rs +++ b/crates/spatialrust-gpu/src/lib.rs @@ -52,9 +52,9 @@ pub use device::{GpuDevice, WgpuDevice}; #[cfg(feature = "gpu-wgpu")] pub use kernels::{ - GpuFpfhDescriptors, GpuFpfhInput, GpuFpfhNeighborhoods, GpuNeighborhoodKind, - GpuNeighborhoodReadback, GpuNeighborhoods, GpuQueryBatch, GPU_FPFH_DESCRIPTOR_LEN, - MAX_GPU_K_NEIGHBORS, + GpuFpfhDescriptors, GpuFpfhInput, GpuFpfhNeighborhoods, GpuIssInput, GpuIssKeypoints, + GpuIssNeighborhoods, 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 317c1d51..14d96186 100644 --- a/crates/spatialrust-gpu/src/runtime.rs +++ b/crates/spatialrust-gpu/src/runtime.rs @@ -26,6 +26,7 @@ pub struct WgpuRuntime { pub(crate) nearest_pipeline: OnceLock, pub(crate) neighborhood_pipelines: OnceLock, pub(crate) fpfh_pipelines: OnceLock, + pub(crate) iss_pipelines: OnceLock, pub(crate) point_estimator_pipelines: OnceLock, max_gather_channels: u32, @@ -390,6 +391,7 @@ impl WgpuRuntime { nearest_pipeline: OnceLock::new(), neighborhood_pipelines: OnceLock::new(), fpfh_pipelines: OnceLock::new(), + iss_pipelines: OnceLock::new(), point_estimator_pipelines: OnceLock::new(), max_gather_channels, upload_pool: GpuBufferPool::default(), diff --git a/crates/spatialrust-gpu/src/shaders/iss.wgsl b/crates/spatialrust-gpu/src/shaders/iss.wgsl new file mode 100644 index 00000000..1e84a487 --- /dev/null +++ b/crates/spatialrust-gpu/src/shaders/iss.wgsl @@ -0,0 +1,167 @@ +struct Params { + count: u32, salient_capacity: u32, non_max_capacity: u32, min_neighbors: u32, + gamma_21: f32, gamma_32: f32, radius: f32, non_max_radius_squared: f32, +} +struct Status { + bad_salient: atomic, bad_non_max: atomic, + required_salient: atomic, required_non_max: atomic, + malformed: atomic, last_bad_plus_one: atomic, + numerical: atomic, reserved: atomic, +} +@group(0) @binding(0) var params: Params; +@group(0) @binding(1) var positions: array>; +@group(0) @binding(2) var salient: array; +@group(0) @binding(3) var non_max: array; +@group(0) @binding(4) var scores: array>; +@group(0) @binding(5) var flags: array; +@group(0) @binding(6) var status: Status; + +// Compensated scalar accumulation reduces covariance summation error. +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 finite(value: f32) -> bool { return abs(value)<=3.4028234663852886e38; } +fn record_bad(i: u32, total: u32, second: bool, malformed: bool) { + if second { + atomicAdd(&status.bad_non_max,1u); atomicMax(&status.required_non_max,total); + } else { + atomicAdd(&status.bad_salient,1u); atomicMax(&status.required_salient,total); + } + if malformed { atomicAdd(&status.malformed,1u); } + atomicMax(&status.last_bad_plus_one,i+1u); +} +fn validate_row(i: u32, second: bool) -> bool { + var capacity=params.salient_capacity; + if second { capacity=params.non_max_capacity; } + let base=i*(2u+2u*capacity); + var stored=0u; var total=0u; + if second { stored=non_max[base]; total=non_max[base+1u]; } + else { stored=salient[base]; total=salient[base+1u]; } + if stored>capacity || totalparams.count { + record_bad(i,total,second,true); return false; + } + if stored!=total { record_bad(i,total,second,false); return false; } + var bound=params.radius*params.radius; + if second { bound=params.non_max_radius_squared; } + for (var n=0u;n(bits); + if j>=params.count || !(distance>=0.0 && distance<=bound) { + record_bad(i,total,second,true); return false; + } + } + return true; +} + +// Max-pivot Jacobi, with the CPU's 32-iteration/physical 1e-12 stopping rule. +// A scaled hypot avoids overflow in the tangent calculation; no eigenvectors +// are computed because ISS only consumes eigenvalues. +fn eigenvalues(matrix: array,3>, epsilon: f32) -> vec3 { + var a=matrix; + 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]; + } + } + } + var values=vec3(a[0][0],a[1][1],a[2][2]); + if values.x>values.y { let t=values.x;values.x=values.y;values.y=t; } + if values.y>values.z { let t=values.y;values.y=values.z;values.z=t; } + if values.x>values.y { let t=values.x;values.x=values.y;values.y=t; } + return values; +} + +@compute @workgroup_size(64) +fn build_saliency(@builtin(global_invocation_id) id: vec3) { + let i=id.x; + if i>=params.count { return; } + scores[i]=vec4(0.0); + // Validate BOTH graphs before any selection pass; no late status races. + let valid_salient=validate_row(i,false); + let valid_non_max=validate_row(i,true); + if !valid_salient || !valid_non_max { return; } + let base=i*(2u+2u*params.salient_capacity); let count=salient[base]; + if count<3u || count,3>; + for (var n=0u;n,6>; + for (var n=0u;n; + for (var axis=0u;axis<3u;axis+=1u) { d[axis]=(relative[axis]-mean[axis].x)-mean[axis].y; } + let products=array(d.x*d.x,d.y*d.y,d.z*d.z,d.x*d.y,d.x*d.z,d.y*d.z); + for (var b=0u;b<6u;b+=1u) { covariance[b]=add(covariance[b],products[b]); } + } + var cov: array; + for (var b=0u;b<6u;b+=1u) { cov[b]=(covariance[b].x+covariance[b].y)/f32(count); } + let matrix=array,3>(vec3(cov[0],cov[3],cov[4]),vec3(cov[3],cov[1],cov[5]),vec3(cov[4],cov[5],cov[2])); + let values=eigenvalues(matrix,1e-12/scale_squared); + let physical=values*scale_squared; + if !finite(physical.x) || !finite(physical.y) || !finite(physical.z) { + atomicAdd(&status.numerical,1u);atomicMax(&status.last_bad_plus_one,i+1u);return; + } + var eligible=0.0; + if values.z>0.0 && values.y>0.0 && values.y/values.z(physical,eligible); +} + +@compute @workgroup_size(64) +fn select_keypoints(@builtin(global_invocation_id) id: vec3) { + let i=id.x; + if i>=params.count { return; } + flags[i]=0u; + if atomicLoad(&status.bad_salient)>0u || atomicLoad(&status.bad_non_max)>0u || atomicLoad(&status.numerical)>0u { + scores[i]=vec4(0.0);return; + } + let score=scores[i]; + if score.w==0.0 { return; } + let base=i*(2u+2u*params.non_max_capacity); let count=non_max[base]; + for (var n=0u;n>()); + assert!(gpu.transfers.host_to_device_bytes() >= 9 * 32); + assert_eq!(gpu.transfers.device_to_host_bytes(), 9 * 4 + 32); + assert!(detector.detect_with_gpu(&runtime, &input, 1, 9).is_err()); +} diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index d41d827f..7e0ea80a 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; ISS/MLS device compute follows. +not change feature algorithm placement; MLS device compute follows. FPFH device computation now consumes this radius-row representation through `GpuFpfhInput` and neighborhoods tied to that exact uploaded input. Explicit @@ -42,6 +42,16 @@ 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. +ISS device computation also consumes input-bound radius rows via `GpuIssInput`. +The uploaded query positions double as covariance samples. Two device passes +validate both graphs, compute compensated local covariance/Jacobi saliency and +select local maxima. Incomplete/malformed graphs or numerical underflow invalidate +all flags/scores. Owned result buffers outlive their inputs; GPU consumers must +inspect status. `feature-iss-gpu` adds explicit detection with transfer receipts, +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. + | Area | Decision | | --- | --- | | Core model | `SpatialTensor` + `SpatialIndex` + `SpatialAlgorithm` + `SpatialRuntime` | diff --git a/docs/FEATURE_MATRIX.md b/docs/FEATURE_MATRIX.md index bfa18541..4278303c 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 +`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 +construction, ordered flag compaction and schema/attribute gathering remain. +Capacities include self and must fit all rows, even rejected candidates. +Positions/extents must be finite, radii need positive finite nonzero squares and +gamma values must be finite in [0,1], stricter than ordinary CPU ISS. Evaluated +nonzero local extents with subnormal/zero f32 squares return numerical errors. +GPU f32 decisions can differ near thresholds/NMS ties from CPU f64. `GpuIssInput` +under `gpu-wgpu` exposes resident graphs/flags and eigen diagnostics. Default +readback is 4 bytes per input point plus 32 status bytes; optional score readback +costs 16 bytes per point plus 32 status bytes. No intermediate graph readback or +automatic CPU fallback occurs. See [ISS evidence](../notes/2026-10-13_gpu_iss_keypoints.md). + `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 @@ -30,7 +44,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 consumes radius rows as described +and preserve query/target semantics. GPU FPFH/ISS 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). @@ -91,6 +105,7 @@ before allocation; no dependency, public signature or image readback is added. | `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 | +| `feature-iss-gpu` | GPU radius graphs, covariance/eigenvalue saliency and local maxima | Explicit detection/receipts; CPU tree, compaction and attribute gathering | | `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 a85e9e3d..2726e2af 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -1128,7 +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 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` | +| 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` | | 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 | +| 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 | | 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_iss_keypoints.md b/notes/2026-10-13_gpu_iss_keypoints.md new file mode 100644 index 00000000..7e011543 --- /dev/null +++ b/notes/2026-10-13_gpu_iss_keypoints.md @@ -0,0 +1,205 @@ +# Resident GPU ISS keypoints + +Date: 2026-10-13. Named slice: ISS radius graphs, covariance/eigenvalue saliency +and local maxima on GPU, following the resident neighborhood and FPFH APIs. + +## Execution and transfer contract + +`feature-iss-gpu` exposes +`IssKeypointDetector::detect_with_gpu(runtime, cloud, salient_capacity, non_max_capacity)` +in features and the facade. Ordinary detection retains its existing CPU path. +The explicit GPU wrapper builds a CPU KD-tree, uploads it and original queries, +runs both radius searches and the two ISS passes on device, then downloads +flags/status. Ordered flag compaction and schema/attribute gathering run on CPU. +No neighborhood/eigenvalue readback or automatic CPU fallback occurs here. + +`GpuIssInput::upload` under `gpu-wgpu` is reusable on one runtime. The query +position buffer is shared with covariance computation; no duplicate sample +buffer is uploaded. 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. Empty +input transfers zero. CPU tree construction is external to this resident API. + +`neighborhoods_resident` retains separate salient and non-max radius graphs tied +to that exact uploaded input. Creating them submits search without readback or +a completion wait. Capacities include self and must fit ALL matches in every +row, even rows that will not qualify as candidates. A truncated prefix is not +accepted as a complete ISS neighborhood. + +`keypoints_resident` uses two ordered device passes in one submission. The first +validates BOTH graphs' counts, target indices and stored distance gates before +any downstream selection. It computes centroid covariance in coordinates +relative to the query, scaled by actual local extent with compensated scalar +accumulation. Max-pivot Jacobi computes ascending eigenvalues without unused +eigenvectors, using 32 iterations and a physical 1e-12 stopping rule. A scaled +hypot avoids tangent overflow. Gamma ratio tests produce eligibility; the second +pass retains local maxima of the smallest eigenvalue. Equal COMPUTED GPU maxima +all survive. Original point ordering is retained in flags and downloaded indices. + +Positions and global extent squares must be finite. Radii must be positive with +finite nonzero f32 squares; gamma values must be finite in [0,1]. These checks +are stricter than ordinary CPU ISS and do not change its configuration behavior. +Minimum support includes self and covariance always needs at least three points. +Huge usize support requirements map to u32::MAX and reject all supported point +counts. Zero local extent produces an ineligible zero row; an evaluated nonzero +extent with a subnormal/zero f32 square sets numerical-error status. A huge search +radius around a small cloud does not dictate covariance scaling. + +Returned `GpuIssKeypoints` owns its flags, scores and status independently of +inputs and other results; none are pooled. Flags are original-order u32 0/1. +Scores are original-order vec4: ascending physical eigenvalues and +eligibility 0/1. Rejected ratios still expose eigenvalues; insufficient support +has a zero score row. Empty physical buffers have unused 4/16 bytes respectively. +The shader requires six storage bindings, one 32-byte uniform and 64 invocations +per workgroup. Buffer sizes, shader indices, dispatch, host bounds and bindings +are checked before allocation. No dependency/core change or unsafe code is added. + +Eight u32 status words contain bad salient rows, bad NMS rows, required salient +capacity, required NMS capacity, malformed rows, last bad row plus one, numerical +errors and a reserved word. Any nonzero status invalidates the complete result. +The second pass zeros ALL flags/scores after graph/numerical errors; downloads +reject status before decoding. Custom GPU consumers must check status and treat +exposed storage as immutable application data. + +Resident search/ISS transfers zero payload bytes. `download` waits and explicitly +copies `4*N + 32` D2H bytes, then compacts flags on CPU. Optional +`download_scores` separately waits and copies `16*N + 32` D2H bytes for diagnostic +eigen/eligibility rows; the wrapper never calls it. Both empty transfers are +zero. Uniform controls are excluded from payload receipts. WASM evidence is +compilation only; browser execution of synchronous downloads was not tested. + +## Precision and correctness evidence + +GPU covariance/eigenvalue arithmetic is f32; CPU ISS uses f64. Thresholds, +nearly zero eigenvalues and nearly tied NMS scores can produce different indices. +This API does not guarantee identical CPU indices across inputs/backends or +equal computed saliencies for all mathematically identical neighborhoods. The +existing CPU implementation remains available when its exact decisions matter. + +Four explicitly run adapter integration tests cover planar, curved/noisy and +coincident-pair fixtures at 0/1/2/17/63/64/65/257/4097 points and radii 0.075/0.16. +They audit 27,396 decisions and 82,188 eigenvalues. Every eigenvalue is compared +with f64 covariance/Jacobi under the fixture bound +`1e-12 + 5e-6 * abs(reference covariance trace)`. Maximum observed absolute +eigenvalue error is 0.000000007803552025051341. This is selected-fixture evidence, +not a universal numerical error guarantee. + +All device flags match independent CPU NMS using the downloaded GPU scores and +CPU radius graph. Compared with ordinary f64 detection, 202 decisions differ +(0.74% of audited input decisions). Every difference is checked against a +threshold uncertainty band or a competing saliency gap within the sum of the +two eigenvalue bounds. Radius-0.16, 4,097-point curved/coincident cases match CPU +indices exactly; radius-0.075 cases and partial narrow clouds include differences. +The accuracy receipts retain every case, including adverse cases. + +Independent analytic box tests verify eigenvalues 0.0004/0.0016/0.0036 within +2e-9 and all eight equal maxima surviving; isolated/linear/coincident support +has no keypoints. Tests exercise both graph overflows and poisoned indices, +exact status and complete output zeroing, empty transfers, numeric underflow, +huge radii, support limits, retained results, four concurrent callers and +runtime/frame ownership. Scalar F32/F64/U8/U16/U32/I32 attributes (including +NaN payloads and signed zero), custom position semantics, schema and metadata +are preserved bit for bit in the selected cloud. + +All-feature GPU regression (203 passed, 23 ignored), feature library regression +(23 passed, two ignored), facade integration, all-feature and CPU-only +Clippy with warnings denied, minimal GPU compilation, WASM feature compilation +and workspace formatting passed. Source comparison verifies that existing CPU +ISS math, selection and gathering are unchanged. Two isolated feature CI cases +are added; 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 omitted. Three warmups per stage precede 20 paired rounds with rotating +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. + +CPU timings use the current public detector, including its bounded worker policy +(four workers at 65,537 points, one for smaller sizes). Full GPU timings include +validation, CPU tree construction, input uploads, BOTH radius graphs, graph +validation, covariance/eigensolves, NMS, flags/status download, completion wait, +ordered compaction and selected-cloud gathering. Runtime/pipeline creation is +warmed outside timers. Result verification/drop is excluded; internal staging +and tree destruction is included in complete calls. + +Resident timings include owned result allocation, graph validation, covariance/ +eigenvalue saliency, NMS, submission and `wait_idle`. They exclude CPU preparation, +input uploads, BOTH neighborhood searches and download/compaction/gathering. +These columns measure a prepared device phase and cannot establish whole-call +or global-registration performance on arbitrary inputs. + +P is a planar 128-column grid with spacing 0.05; C adds curved height and +deterministic z noise; D duplicates consecutive C points. Salient radius is +0.11 for P and 0.16 for C/D; NMS radius is 0.75 times that. Salient capacity is +64 for P/C and 96 for D; NMS capacity is 64 throughout. Numerical audits run +before timing; repeated CPU calls match their CPU reference, and repeated GPU +calls match their separately audited GPU reference. Selected counts and +CPU/GPU decision differences accompany performance: the two precision modes +can return different indices, so ratios 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 with seeds 20261013 plus zero-based record index. +They do not establish a deployment crossover or remove run-to-run variability. + +| Survey | Points | Fixture | Workers | CPU ms | Full GPU ms | Resident ms | CPU/GPU selected | Decision differences | Paired full-call ratio (95% interval) | +| --- | ---: | :---: | ---: | ---: | ---: | ---: | ---: | ---: | --- | +| 1 | 257 | P | 1 | 0.250 | 1.465 | 0.201 | 256/256 | 0 | 0.17x (0.16-0.18) | +| 1 | 257 | C | 1 | 0.389 | 1.492 | 0.242 | 44/30 | 14 | 0.28x (0.26-0.30) | +| 1 | 257 | D | 1 | 0.322 | 1.436 | 0.200 | 248/248 | 0 | 0.23x (0.22-0.24) | +| 1 | 4,097 | P | 1 | 5.344 | 3.111 | 0.241 | 624/624 | 0 | 1.72x (1.66-1.80) | +| 1 | 4,097 | C | 1 | 12.167 | 3.006 | 0.298 | 421/421 | 0 | 3.99x (3.81-4.18) | +| 1 | 4,097 | D | 1 | 14.596 | 2.969 | 0.395 | 444/444 | 0 | 4.96x (4.77-5.16) | +| 1 | 16,384 | P | 1 | 23.620 | 7.394 | 0.548 | 1008/1008 | 0 | 3.13x (2.71-3.56) | +| 1 | 16,384 | C | 1 | 57.514 | 24.491 | 10.328 | 2046/2046 | 0 | 2.74x (2.50-3.01) | +| 1 | 16,384 | D | 1 | 87.431 | 23.404 | 10.637 | 1762/1762 | 0 | 3.06x (2.67-3.48) | +| 1 | 65,537 | P | 4 | 39.465 | 36.799 | 9.382 | 2544/2544 | 0 | 1.06x (1.00-1.12) | +| 1 | 65,537 | C | 4 | 97.191 | 32.245 | 3.352 | 5412/5412 | 0 | 2.69x (2.37-3.05) | +| 1 | 65,537 | D | 4 | 125.311 | 35.082 | 8.125 | 6483/6483 | 0 | 3.55x (3.32-3.76) | +| 2 | 257 | P | 1 | 0.221 | 1.296 | 0.165 | 256/256 | 0 | 0.17x (0.16-0.18) | +| 2 | 257 | C | 1 | 0.341 | 1.193 | 0.167 | 44/30 | 14 | 0.28x (0.27-0.29) | +| 2 | 257 | D | 1 | 0.273 | 1.213 | 0.176 | 248/248 | 0 | 0.22x (0.21-0.23) | +| 2 | 4,097 | P | 1 | 4.002 | 2.484 | 0.241 | 624/624 | 0 | 1.68x (1.64-1.72) | +| 2 | 4,097 | C | 1 | 13.223 | 3.265 | 0.323 | 421/421 | 0 | 3.82x (3.66-3.98) | +| 2 | 4,097 | D | 1 | 12.090 | 3.077 | 0.530 | 444/444 | 0 | 2.87x (2.06-3.97) | +| 2 | 16,384 | P | 1 | 19.347 | 6.616 | 0.525 | 1008/1008 | 0 | 2.91x (2.79-3.03) | +| 2 | 16,384 | C | 1 | 68.500 | 14.967 | 4.334 | 2046/2046 | 0 | 5.65x (4.82-6.66) | +| 2 | 16,384 | D | 1 | 82.495 | 22.265 | 9.680 | 1762/1762 | 0 | 3.12x (2.75-3.51) | +| 2 | 65,537 | P | 4 | 43.909 | 40.607 | 9.406 | 2544/2544 | 0 | 1.06x (0.97-1.14) | +| 2 | 65,537 | C | 4 | 95.802 | 32.239 | 3.344 | 5412/5412 | 0 | 2.68x (2.33-3.07) | +| 2 | 65,537 | D | 4 | 130.787 | 37.264 | 8.104 | 6483/6483 | 0 | 3.59x (3.40-3.79) | + +At 65,537 points, both surveys have identical CPU/GPU indices on C/D, with +paired complete-call ratios 2.68-3.59x against four-worker CPU detection. P is +about 1.06x and its bootstrap intervals reach/include 1, giving no clear gain. +All 257-point cases favor CPU, and C differs by 14 of 257 decisions. Timing +variation remains visible between surveys and stages; these synthetic fixtures +establish neither a universal gain nor an accuracy-equivalent deployment claim. +Benchmark audits have maximum absolute eigenvalue error 0.00000000935703506517449 +under the same trace-based fixture bound; diagnostics run outside timed calls. + +## Reproduction + +From `C:/Users/rsasa/Workspace/SpatialRust`, run adapter checks explicitly: + +```powershell +cargo test -p spatialrust-features --no-default-features --features feature-iss-gpu --test gpu_iss -- --include-ignored --nocapture +cargo test -p spatialrust --no-default-features --features feature-iss-gpu --test gpu_iss_smoke -- --include-ignored +``` + +Run this exact benchmark filter alone twice to exclude the included reference +module's adapter tests from timing runs: + +```powershell +cargo test -p spatialrust-features --no-default-features --features feature-iss-gpu --test gpu_iss_benchmark --release completed_iss_calls -- --include-ignored --nocapture +``` + +Local untracked evidence: + +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-iss-integration-tests.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-iss-native-tests.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-iss-feature-tests.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-iss-facade-tests.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-iss-benchmark-1.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-iss-benchmark-2.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-iss-measurements.json`