diff --git a/CHANGELOG.md b/CHANGELOG.md index 839e6640..d82c0470 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,12 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`). ### Changed +- ISS keypoint detection shares one KD-tree across bounded scoped CPU workers + for covariance saliency and local-maximum selection under `feature-iss`. + Saliency completes before selection; small/one-worker clouds stay serial. + Equal maxima, selected point order, schema, attributes and metadata remain + unchanged. Parallel selection uses one temporary byte per point, without a + new dependency or public signature. - Point-cloud boundary detection shares one KD-tree across bounded scoped CPU workers under `feature-boundary`. Workers reuse angular scratch and write disjoint mask bands directly; small/one-worker clouds stay serial. Boundary diff --git a/crates/spatialrust-features/src/iss.rs b/crates/spatialrust-features/src/iss.rs index de3bdecc..0d9a7dd7 100644 --- a/crates/spatialrust-features/src/iss.rs +++ b/crates/spatialrust-features/src/iss.rs @@ -10,7 +10,7 @@ use spatialrust_core::{ HasPositions3, PointBuffer, PointBufferSet, PointCloud, SpatialError, SpatialResult, }; use spatialrust_math::{symmetric_eigen3, Mat3}; -use spatialrust_search::{KdTree, RadiusSearchIndex}; +use spatialrust_search::{parallel_worker_count, KdTree, RadiusSearchIndex}; /// Configuration for [`IssKeypointDetector`]. #[derive(Clone, Copy, Debug, PartialEq)] @@ -76,7 +76,19 @@ impl IssKeypointDetector { } /// Detects keypoints and returns their indices and a keypoint sub-cloud. + /// + /// Shares one KD-tree across bounded scoped CPU workers. Saliency completes + /// before local-maximum selection reads it. Small/one-worker calls run on + /// the caller; neighborhood arithmetic, ties and point order are preserved. pub fn detect(&self, input: &PointCloud) -> SpatialResult { + self.detect_with_workers(input, parallel_worker_count(input.len())) + } + + fn detect_with_workers( + &self, + input: &PointCloud, + workers: usize, + ) -> SpatialResult { if self.config.salient_radius <= 0.0 || self.config.non_max_radius <= 0.0 { return Err(SpatialError::InvalidArgument("ISS radii must be positive".to_owned())); } @@ -84,45 +96,67 @@ impl IssKeypointDetector { let (x, y, z) = input.positions3()?; let len = input.len(); let tree = KdTree::from_slices(x, y, z); + let workers = workers.max(1).min(len.max(1)); // Saliency = smallest eigenvalue (lambda3) where the eigenvalue-ratio // tests pass; NaN marks a non-salient point. let mut saliency = vec![f32::NAN; len]; - for i in 0..len { - let neighbors = tree.radius_search(x[i], y[i], z[i], self.config.salient_radius); - if neighbors.len() < self.config.min_neighbors { - continue; - } - let Some(eigenvalues) = neighborhood_eigenvalues(x, y, z, &neighbors) else { - continue; - }; - // Ascending order: l3 <= l2 <= l1. - let (l3, l2, l1) = (eigenvalues[0], eigenvalues[1], eigenvalues[2]); - if l1 <= 0.0 || l2 <= 0.0 { - continue; - } - if (l2 / l1) < f64::from(self.config.gamma_21) - && (l3 / l2) < f64::from(self.config.gamma_32) - { - saliency[i] = l3 as f32; + fill_bands(&mut saliency, workers, |start, band| { + for (offset, value) in band.iter_mut().enumerate() { + let i = start + offset; + let neighbors = tree.radius_search(x[i], y[i], z[i], self.config.salient_radius); + if neighbors.len() < self.config.min_neighbors { + continue; + } + let Some(eigenvalues) = neighborhood_eigenvalues(x, y, z, &neighbors) else { + continue; + }; + // Ascending order: l3 <= l2 <= l1. + let (l3, l2, l1) = (eigenvalues[0], eigenvalues[1], eigenvalues[2]); + if l1 <= 0.0 || l2 <= 0.0 { + continue; + } + if (l2 / l1) < f64::from(self.config.gamma_21) + && (l3 / l2) < f64::from(self.config.gamma_32) + { + *value = l3 as f32; + } } - } + }); - // Non-maximum suppression: keep points whose saliency is the strict - // maximum within `non_max_radius`. - let mut indices = Vec::new(); - for i in 0..len { + // The completed saliency array is immutable throughout this phase. + // Equal maxima all survive, retaining the original >= tie behavior. + let is_keypoint = |i: usize| { let s = saliency[i]; if s.is_nan() { - continue; + return false; } let neighbors = tree.radius_search(x[i], y[i], z[i], self.config.non_max_radius); - let is_local_max = neighbors.iter().all(|n| { + neighbors.iter().all(|n| { let other = saliency[n.index]; n.index == i || other.is_nan() || s >= other + }) + }; + let mut indices = Vec::new(); + if workers == 1 { + for i in 0..len { + if is_keypoint(i) { + indices.push(i); + } + } + } else { + // One temporary byte per point avoids worker-local index vectors + // and a second result-sized index allocation during ordered joins. + let mut selected = vec![false; len]; + fill_bands(&mut selected, workers, |start, band| { + for (offset, value) in band.iter_mut().enumerate() { + *value = is_keypoint(start + offset); + } }); - if is_local_max { - indices.push(i); + for (i, selected) in selected.into_iter().enumerate() { + if selected { + indices.push(i); + } } } @@ -131,6 +165,34 @@ impl IssKeypointDetector { } } +/// Fills disjoint output bands, joining all workers before the next phase. +fn fill_bands(output: &mut [T], workers: usize, fill: F) { + if workers == 1 { + fill(0, output); + return; + } + let chunk_size = output.len().div_ceil(workers); + std::thread::scope(|scope| { + for (band, chunk) in output.chunks_mut(chunk_size).enumerate() { + let fill = &fill; + scope.spawn(move || fill(band * chunk_size, chunk)); + } + }); +} + +#[cfg(all(test, not(target_arch = "wasm32")))] +#[path = "../tests/fixtures/iss_serial.rs"] +#[allow(dead_code)] +mod previous; + +#[cfg(all(test, not(target_arch = "wasm32")))] +#[path = "iss_parallel_tests.rs"] +mod parallel_tests; + +#[cfg(all(test, not(target_arch = "wasm32")))] +#[path = "iss_parallel_benchmark.rs"] +mod parallel_benchmark; + /// Eigenvalues (ascending) of the neighborhood covariance about its centroid. fn neighborhood_eigenvalues( x: &[f32], diff --git a/crates/spatialrust-features/src/iss_parallel_benchmark.rs b/crates/spatialrust-features/src/iss_parallel_benchmark.rs new file mode 100644 index 00000000..efeb432c --- /dev/null +++ b/crates/spatialrust-features/src/iss_parallel_benchmark.rs @@ -0,0 +1,54 @@ +use super::{ + parallel_tests::{assert_result, fixture, previous_result}, + *, +}; +use std::{hint::black_box, time::Instant}; + +#[test] +#[ignore = "completed ISS detection comparison; explicit native invocation only"] +fn paired_previous_serial_parallel_iss() { + let available = std::thread::available_parallelism().map_or(1, |n| n.get()); + eprintln!("available_parallelism={available}; 2 warmups + 20 rotating rounds; default shared staging policy"); + eprintln!("Timers include validation, KD-tree build, saliency allocation/covariance, local-maximum selection, worker spawn/join, ordered indices, schema/attribute gathering and scratch/tree destruction; exclude fixture creation, result verification and returned result drop."); + eprintln!("points,kind,radius,workers,round,previous_ms,serial_ms,parallel_ms"); + for points in [512usize, 16384, 32768, 65536] { + for kind in 0..3 { + let radius = if kind == 0 { 0.11 } else { 0.16 }; + let config = IssKeypointConfig::with_radii(radius, radius * 0.75); + let cloud = fixture(points, kind); + let expected = previous_result(&cloud, config); + let detector = IssKeypointDetector::new(config); + let workers = parallel_worker_count(points); + let run = |mode| { + if mode == 0 { + let previous = previous_result(&cloud, config); + IssKeypointResult { indices: previous.indices, keypoints: previous.keypoints } + } else if mode == 1 { + detector.detect_with_workers(&cloud, 1).unwrap() + } else { + detector.detect(&cloud).unwrap() + } + }; + for _ in 0..2 { + for mode in 0..3 { + assert_result(&run(mode), &expected); + } + } + for round in 0..20 { + let mut ms = [0.0; 3]; + for offset in 0..3 { + let mode = (round + offset) % 3; + let start = Instant::now(); + let actual = run(mode); + ms[mode] = start.elapsed().as_secs_f64() * 1000.0; + assert_result(&actual, &expected); + black_box(&actual); + } + eprintln!( + "{points},{kind},{radius},{workers},{round},{:.6},{:.6},{:.6}", + ms[0], ms[1], ms[2] + ); + } + } + } +} diff --git a/crates/spatialrust-features/src/iss_parallel_tests.rs b/crates/spatialrust-features/src/iss_parallel_tests.rs new file mode 100644 index 00000000..213db3d5 --- /dev/null +++ b/crates/spatialrust-features/src/iss_parallel_tests.rs @@ -0,0 +1,250 @@ +use super::*; +use spatialrust_core::{ + DType, FieldSemantic, PointCloudBuilder, PointField, PointSchema, SpatialMetadata, + StandardSchemas, Timestamp, +}; + +pub(super) fn fixture(len: usize, kind: usize) -> PointCloud { + let mut columns: [Vec; 3] = std::array::from_fn(|_| Vec::with_capacity(len)); + for i in 0..len { + 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); + } + } + let mut buffers = PointBufferSet::new(); + for (name, column) in ["x", "y", "z"].into_iter().zip(columns) { + buffers.insert(name, PointBuffer::F32(column)); + } + PointCloud::try_from_parts(StandardSchemas::point_xyz(), buffers, Default::default()).unwrap() +} + +fn previous_detector(config: IssKeypointConfig) -> previous::IssKeypointDetector { + previous::IssKeypointDetector::new(previous::IssKeypointConfig { + salient_radius: config.salient_radius, + non_max_radius: config.non_max_radius, + gamma_21: config.gamma_21, + gamma_32: config.gamma_32, + min_neighbors: config.min_neighbors, + }) +} + +pub(super) fn previous_result( + cloud: &PointCloud, + config: IssKeypointConfig, +) -> previous::IssKeypointResult { + previous_detector(config).detect(cloud).unwrap() +} + +pub(super) fn assert_result(actual: &IssKeypointResult, expected: &previous::IssKeypointResult) { + assert_eq!(actual.indices, expected.indices); + assert_eq!(actual.keypoints, expected.keypoints); +} + +#[test] +fn frozen_selection_and_clouds_match_across_worker_bands_and_public_staging() { + let mut decisions = 0; + for len in [0usize, 1, 17, 4095, 4096, 4103, 32771] { + for kind in 0..3 { + let cloud = fixture(len, kind); + for radius in [0.075, 0.16] { + let config = IssKeypointConfig { + min_neighbors: if kind == 2 { 0 } else { 5 }, + ..IssKeypointConfig::with_radii(radius, radius * 0.75) + }; + let expected = previous_result(&cloud, config); + for workers in [1, 2, 6] { + let actual = IssKeypointDetector::new(config) + .detect_with_workers(&cloud, workers) + .unwrap(); + assert_result(&actual, &expected); + decisions += len; + } + assert_result(&IssKeypointDetector::new(config).detect(&cloud).unwrap(), &expected); + decisions += len; + } + } + } + eprintln!("Verified {decisions} exact ISS selection decisions and selected clouds against frozen serial detection."); +} + +#[test] +fn independent_covariance_ties_and_degenerate_support_are_preserved() { + let mut builder = PointCloudBuilder::xyz(); + for x in [-0.06, 0.06] { + for y in [-0.04, 0.04] { + for z in [-0.02, 0.02] { + builder.push_point([x, y, z]).unwrap(); + } + } + } + // An isolated point does not have enough covariance support. + builder.push_point([100.0, 100.0, 100.0]).unwrap(); + let cloud = builder.build().unwrap(); + let (x, y, z) = cloud.positions3().unwrap(); + let tree = KdTree::from_slices(x, y, z); + let eigenvalues = + neighborhood_eigenvalues(x, y, z, &tree.radius_search(x[0], y[0], z[0], 0.2)).unwrap(); + for (actual, extent) in eigenvalues.into_iter().zip([0.02_f32, 0.04, 0.06]) { + assert!((actual - f64::from(extent).powi(2)).abs() < 1e-14); + } + let config = IssKeypointConfig::with_radii(0.2, 1.0); + for workers in [0, 1, 2, 6, usize::MAX] { + let actual = IssKeypointDetector::new(config).detect_with_workers(&cloud, workers).unwrap(); + assert_eq!(actual.indices, (0..8).collect::>(), "equal maxima must all survive"); + assert_eq!(actual.keypoints.len(), 8); + } + for points in [ + vec![[0.0, 0.0, 0.0]; 33], + (0..33).map(|i| [i as f32 * 0.01, 0.0, 0.0]).collect(), + vec![[-0.0, 0.0, -0.0], [100.0, 0.0, 0.0]], + ] { + let mut builder = PointCloudBuilder::xyz(); + for p in points { + builder.push_point(p).unwrap(); + } + let cloud = builder.build().unwrap(); + for min_neighbors in [0, 3, usize::MAX] { + let config = + IssKeypointConfig { min_neighbors, ..IssKeypointConfig::with_radii(0.025, 0.02) }; + for workers in [1, 2, 6] { + let actual = + IssKeypointDetector::new(config).detect_with_workers(&cloud, workers).unwrap(); + assert!(actual.indices.is_empty()); + assert_result(&actual, &previous_result(&cloud, config)); + } + } + } +} + +#[test] +fn selected_indices_preserve_schema_metadata_and_all_attribute_words() { + let source = fixture(32771, 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 config = IssKeypointConfig::with_radii(0.16, 0.12); + let expected = previous_result(&input, config); + assert!(!expected.indices.is_empty()); + let detector = IssKeypointDetector::new(config); + for result in + [detector.detect_with_workers(&input, 6).unwrap(), detector.detect(&input).unwrap()] + { + assert_eq!(result.indices, expected.indices); + assert!(result.indices.windows(2).all(|w| w[0] < w[1])); + assert_eq!(result.keypoints.schema(), input.schema()); + assert_eq!(result.keypoints.metadata(), input.metadata()); + for field in input.schema().fields() { + match ( + result.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] +fn validation_unusual_configuration_and_concurrent_calls_match_serial_contract() { + let cloud = fixture(17, 1); + let empty = fixture(0, 0); + let singleton = fixture(1, 0); + for input in [&cloud, &empty, &singleton] { + for config in [ + IssKeypointConfig::with_radii(0.0, 0.1), + IssKeypointConfig::with_radii(0.1, -1.0), + IssKeypointConfig::with_radii(f32::NAN, 0.1), + IssKeypointConfig::with_radii(0.1, f32::NAN), + IssKeypointConfig::with_radii(f32::INFINITY, f32::INFINITY), + IssKeypointConfig { gamma_21: f32::NAN, ..IssKeypointConfig::with_radii(0.16, 0.12) }, + IssKeypointConfig { gamma_32: -1.0, ..IssKeypointConfig::with_radii(0.16, 0.12) }, + ] { + let expected = previous_detector(config) + .detect(input) + .map(|r| (r.indices, r.keypoints)) + .map_err(|e| e.to_string()); + for workers in [0, 1, 2, 6, usize::MAX] { + let actual = IssKeypointDetector::new(config) + .detect_with_workers(input, workers) + .map(|r| (r.indices, r.keypoints)) + .map_err(|e| e.to_string()); + assert_eq!(actual, expected); + } + } + } + let config = IssKeypointConfig::with_radii(0.16, 0.12); + let cloud = fixture(4103, 1); + let expected = previous_result(&cloud, config); + std::thread::scope(|scope| { + let calls: Vec<_> = [2, 6] + .into_iter() + .map(|workers| { + let cloud = &cloud; + scope.spawn(move || { + IssKeypointDetector::new(config).detect_with_workers(cloud, workers).unwrap() + }) + }) + .collect(); + for call in calls { + assert_result(&call.join().unwrap(), &expected); + } + }); +} diff --git a/crates/spatialrust-features/tests/fixtures/iss_serial.rs b/crates/spatialrust-features/tests/fixtures/iss_serial.rs new file mode 100644 index 00000000..063ab3ca --- /dev/null +++ b/crates/spatialrust-features/tests/fixtures/iss_serial.rs @@ -0,0 +1,196 @@ +//! Intrinsic Shape Signatures (ISS) keypoint detection. +//! +//! ISS picks a sparse set of geometrically salient points — corners and other +//! spots with variation in all three directions — by thresholding the ratios of +//! the eigenvalues of each point's neighborhood covariance, then keeping local +//! maxima of the smallest eigenvalue. The keypoints are a natural front-end for +//! descriptor-based global registration: far fewer points to describe and match. + +use spatialrust_core::{ + HasPositions3, PointBuffer, PointBufferSet, PointCloud, SpatialError, SpatialResult, +}; +use spatialrust_math::{symmetric_eigen3, Mat3}; +use spatialrust_search::{KdTree, RadiusSearchIndex}; + +/// Configuration for [`IssKeypointDetector`]. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct IssKeypointConfig { + /// Radius of the neighborhood used to build each covariance matrix. + pub salient_radius: f32, + /// Radius for non-maximum suppression of the saliency (smallest eigenvalue). + pub non_max_radius: f32, + /// Upper bound on `lambda2 / lambda1` (rejects flat / planar points). + pub gamma_21: f32, + /// Upper bound on `lambda3 / lambda2` (rejects edge / linear points). + pub gamma_32: f32, + /// Minimum neighbors within `salient_radius` for a point to be considered. + pub min_neighbors: usize, +} + +impl Default for IssKeypointConfig { + fn default() -> Self { + Self { + salient_radius: 0.2, + non_max_radius: 0.15, + gamma_21: 0.975, + gamma_32: 0.975, + min_neighbors: 5, + } + } +} + +impl IssKeypointConfig { + /// Creates a config from the salient and non-max-suppression radii. + #[must_use] + pub fn with_radii(salient_radius: f32, non_max_radius: f32) -> Self { + Self { salient_radius, non_max_radius, ..Self::default() } + } +} + +/// Result of ISS keypoint detection. +#[derive(Clone, Debug, PartialEq)] +pub struct IssKeypointResult { + /// Indices of the keypoints in the input cloud. + pub indices: Vec, + /// Sub-cloud containing only the keypoints (input schema preserved). + pub keypoints: PointCloud, +} + +/// Intrinsic Shape Signatures keypoint detector. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct IssKeypointDetector { + config: IssKeypointConfig, +} + +impl IssKeypointDetector { + /// Creates a detector from config. + #[must_use] + pub const fn new(config: IssKeypointConfig) -> Self { + Self { config } + } + + /// Returns the detector config. + #[must_use] + pub const fn config(&self) -> IssKeypointConfig { + self.config + } + + /// Detects keypoints and returns their indices and a keypoint sub-cloud. + pub fn detect(&self, input: &PointCloud) -> SpatialResult { + if self.config.salient_radius <= 0.0 || self.config.non_max_radius <= 0.0 { + return Err(SpatialError::InvalidArgument("ISS radii must be positive".to_owned())); + } + + let (x, y, z) = input.positions3()?; + let len = input.len(); + let tree = KdTree::from_slices(x, y, z); + + // Saliency = smallest eigenvalue (lambda3) where the eigenvalue-ratio + // tests pass; NaN marks a non-salient point. + let mut saliency = vec![f32::NAN; len]; + for i in 0..len { + let neighbors = tree.radius_search(x[i], y[i], z[i], self.config.salient_radius); + if neighbors.len() < self.config.min_neighbors { + continue; + } + let Some(eigenvalues) = neighborhood_eigenvalues(x, y, z, &neighbors) else { + continue; + }; + // Ascending order: l3 <= l2 <= l1. + let (l3, l2, l1) = (eigenvalues[0], eigenvalues[1], eigenvalues[2]); + if l1 <= 0.0 || l2 <= 0.0 { + continue; + } + if (l2 / l1) < f64::from(self.config.gamma_21) + && (l3 / l2) < f64::from(self.config.gamma_32) + { + saliency[i] = l3 as f32; + } + } + + // Non-maximum suppression: keep points whose saliency is the strict + // maximum within `non_max_radius`. + let mut indices = Vec::new(); + for i in 0..len { + let s = saliency[i]; + if s.is_nan() { + continue; + } + let neighbors = tree.radius_search(x[i], y[i], z[i], self.config.non_max_radius); + let is_local_max = neighbors.iter().all(|n| { + let other = saliency[n.index]; + n.index == i || other.is_nan() || s >= other + }); + if is_local_max { + indices.push(i); + } + } + + let keypoints = gather_indices(input, &indices)?; + Ok(IssKeypointResult { indices, keypoints }) + } +} + +/// Eigenvalues (ascending) of the neighborhood covariance about its centroid. +fn neighborhood_eigenvalues( + x: &[f32], + y: &[f32], + z: &[f32], + neighbors: &[spatialrust_search::Neighbor], +) -> Option<[f64; 3]> { + let count = neighbors.len() as f64; + if count < 3.0 { + return None; + } + let (mut mx, mut my, mut mz) = (0.0_f64, 0.0_f64, 0.0_f64); + for n in neighbors { + mx += f64::from(x[n.index]); + my += f64::from(y[n.index]); + mz += f64::from(z[n.index]); + } + mx /= count; + my /= count; + mz /= count; + + let (mut c00, mut c11, mut c22) = (0.0_f64, 0.0, 0.0); + let (mut c01, mut c02, mut c12) = (0.0_f64, 0.0, 0.0); + for n in neighbors { + let dx = f64::from(x[n.index]) - mx; + let dy = f64::from(y[n.index]) - my; + let dz = f64::from(z[n.index]) - mz; + c00 += dx * dx; + c11 += dy * dy; + c22 += dz * dz; + c01 += dx * dy; + c02 += dx * dz; + c12 += dy * dz; + } + let inv = 1.0 / count; + let covariance = Mat3::::from_rows( + [c00 * inv, c01 * inv, c02 * inv], + [c01 * inv, c11 * inv, c12 * inv], + [c02 * inv, c12 * inv, c22 * inv], + ); + Some(symmetric_eigen3(covariance).eigenvalues) +} + +/// Gathers the selected indices into a new cloud, preserving schema. +fn gather_indices(input: &PointCloud, indices: &[usize]) -> SpatialResult { + let mut buffers = PointBufferSet::new(); + for field in input.schema().fields() { + let source = input.field(&field.name)?; + buffers.insert(field.name.clone(), gather_buffer(source, indices)); + } + PointCloud::try_from_parts(input.schema().clone(), buffers, input.metadata().clone()) +} + +fn gather_buffer(source: &PointBuffer, indices: &[usize]) -> PointBuffer { + match source { + PointBuffer::F32(v) => PointBuffer::from_f32(indices.iter().map(|&i| v[i]).collect()), + PointBuffer::F64(v) => PointBuffer::F64(indices.iter().map(|&i| v[i]).collect()), + PointBuffer::U8(v) => PointBuffer::U8(indices.iter().map(|&i| v[i]).collect()), + PointBuffer::U16(v) => PointBuffer::U16(indices.iter().map(|&i| v[i]).collect()), + PointBuffer::U32(v) => PointBuffer::U32(indices.iter().map(|&i| v[i]).collect()), + PointBuffer::I32(v) => PointBuffer::I32(indices.iter().map(|&i| v[i]).collect()), + } +} diff --git a/crates/spatialrust/tests/iss_parallel_smoke.rs b/crates/spatialrust/tests/iss_parallel_smoke.rs new file mode 100644 index 00000000..b55a40b8 --- /dev/null +++ b/crates/spatialrust/tests/iss_parallel_smoke.rs @@ -0,0 +1,23 @@ +#![cfg(feature = "feature-iss")] + +use spatialrust::{IssKeypointConfig, IssKeypointDetector, PointCloudBuilder}; + +#[test] +fn facade_iss_preserves_equal_maxima_order_and_selected_cloud() { + let mut builder = PointCloudBuilder::xyz(); + for x in [-0.06, 0.06] { + for y in [-0.04, 0.04] { + for z in [-0.02, 0.02] { + builder.push_point([x, y, z]).unwrap(); + } + } + } + builder.push_point([100.0, 100.0, 100.0]).unwrap(); + let input = builder.build().unwrap(); + let detector = IssKeypointDetector::new(IssKeypointConfig::with_radii(0.2, 1.0)); + let result = detector.detect(&input).unwrap(); + assert_eq!(result.indices, (0..8).collect::>()); + assert_eq!(result.keypoints.len(), 8); + assert_eq!(result.keypoints.schema(), input.schema()); + assert_eq!(result.keypoints.metadata(), input.metadata()); +} diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index 57a43f91..6e28d66b 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -313,6 +313,17 @@ selected points in original index order and preserves all field types and metadata. Small/one-worker calls remain serial. Neighbor/angle scratch scales with radius support per active worker; the worker cap is not a scratch-byte or process-wide CPU budget. No device copy or new dependency is introduced. +ISS keypoint detection under `spatialrust-features/feature-iss` uses the same +bounded staging policy in two scoped phases sharing one KD-tree. Workers fill +disjoint slices of the existing four-byte-per-point saliency array; the first +scope joins before local-maximum selection can read any row. Parallel selection +writes a temporary one-byte-per-point mask and scans it in original index order, +dropping it before schema/attribute gathering. Serial calls retain direct index +collection and allocate no selection mask. Unsorted neighborhood order, f64 +centroid/covariance arithmetic, eigenvalue ratios and the existing >= tie rule +are preserved. Per-query neighbor vectors scale with radius support; worker +caps remain per call rather than a scratch-byte or process-wide CPU budget. +No floating-point reduction spans workers and no device transfer is introduced. `spatialrust-image-io` depends on storage, never the reverse; standard codecs are additive, while TIFF and OpenEXR remain independently gated. Calibration datasets, robust solver controls, and residual receipts live in diff --git a/docs/FEATURE_MATRIX.md b/docs/FEATURE_MATRIX.md index 4fd7578a..797bde79 100644 --- a/docs/FEATURE_MATRIX.md +++ b/docs/FEATURE_MATRIX.md @@ -212,6 +212,12 @@ workspace because its build requires a Python toolchain. decisions, selected attributes and metadata are preserved. Small/one-worker calls stay serial, with no additional dependency or device transfer. See [boundary execution measurements](../notes/2026-10-12_parallel_boundary_detection.md). +- `feature-iss` shares one KD-tree across bounded covariance and local-maximum + phases. Saliency joins before selection reads it; parallel selection uses one + temporary byte per point. Small/one-worker calls keep direct serial index + collection. Equal maxima, selected point order, schema, all attribute types + and metadata are preserved, without another dependency or device transfer. + See [ISS detection measurements](../notes/2026-10-12_parallel_iss_keypoints.md). - `GpuStageProfiler` (under `gpu-wgpu`) records named pass timestamps through `run_gpu_vision_chain_profiled` (`gpu-image`) and `run_aoso_voxel_normal_frame_profiled` (`gpu-aoso-staging`). Create an optional diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index a794dfab..dd6d8844 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -1096,6 +1096,7 @@ v1.3 admission gates. | Bounded CPU MLS smoothing | Complete | Independent point fits on shared KD-tree with existing bounded staging and disjoint ordered output bands | Frozen exact-coordinate parity with 1/2/6 workers, curved/duplicate/degenerate/support cases, all attribute types and metadata, feature/facade/WASM checks and twelve completed-position-API comparisons; `notes/2026-10-12_parallel_mls_smoothing.md` | | Bounded CPU FPFH descriptors | Complete | Shared KD-tree, scoped SPFH/weighting phases and direct disjoint descriptor writes using existing point staging | 18,402,648 frozen exact descriptor words with 1/2/6 workers and public staging, independent planar/degenerate cases, validation/concurrency/seeded alignment, feature/facade/WASM checks and twelve completed-descriptor-API comparisons; `notes/2026-10-12_parallel_fpfh_descriptors.md` | | Bounded CPU boundary detection | Complete | Shared KD-tree, reusable worker angle scratch and direct disjoint mask writes using existing point staging | 1,081,992 frozen exact decisions with 1/2/6 workers and public staging, independent complete rim/degenerate/normal cases, all attribute types and metadata, validation/concurrency, feature/facade/WASM checks and twelve completed-mask-API comparisons; `notes/2026-10-12_parallel_boundary_detection.md` | +| Bounded CPU ISS keypoints | Complete | Shared KD-tree, scoped covariance/selection phases and ordered selection with existing point staging | 1,081,992 frozen exact selection decisions and selected clouds with 1/2/6 workers and public staging, independent covariance/tied maxima/degenerate support, all attribute types and metadata, validation/concurrency, feature/facade/WASM checks and twelve completed-detection-API comparisons; `notes/2026-10-12_parallel_iss_keypoints.md` | | GPU spatial execution guards | Complete | Pre-allocation runtime/dimension/device/profiler checks for nearest resize and Sobel, including checked resize coordinate products | CPU/frozen full-RGBA parity, invalid-limit/capacity and ROI/concurrent/live-output/pool checks, image/point regressions and feature/WASM builds; `notes/2026-10-12_gpu_spatial_dispatch_validation.md` | | Bounded GPU gray dispatch | Complete | Balanced multi-row flat workgroups beyond the single-axis limit, independently lazy kernels and pre-allocation arithmetic/device/profiler validation | All 16,777,216 RGB values, 120 forced grid mappings, maximum-counter device probes, ownership/capacity/ROI/concurrency/reuse checks, native regressions and resident timing; `notes/2026-10-12_bounded_gpu_gray_dispatch.md` | | GPU AI execution guards | Complete | Pre-allocation image/index/buffer/dispatch and profiler validation while preserving the original CHW shader | 23,857,560 frozen output words across 192 cases, ROI/concurrency/live tensor reuse, negative-limit/capacity checks, 4096x4096 six-stage chain and feature/native/WASM validation; `notes/2026-10-12_validated_gpu_ai_packing.md` | diff --git a/notes/2026-10-12_parallel_iss_keypoints.md b/notes/2026-10-12_parallel_iss_keypoints.md new file mode 100644 index 00000000..393a355b --- /dev/null +++ b/notes/2026-10-12_parallel_iss_keypoints.md @@ -0,0 +1,157 @@ +# Bounded CPU ISS keypoint detection + +Date: 2026-10-12. This named slice follows +[PR #148](https://github.com/rsasaki0109/SpatialRust/pull/148), source +`e5b0c6a6299350679fc40d521844561c41510b26`. + +## Result and execution boundary + +ISS previously evaluated every neighborhood covariance and local maximum on +the caller. Under `feature-iss`, `IssKeypointDetector::detect` now shares one +caller-built KD-tree across bounded scoped CPU workers. The first phase fills +disjoint contiguous slices of the existing saliency array. Its scope joins +completely before local-maximum selection reads any saliency value. Parallel +selection fills a temporary mask, then gathers indices in ascending original +order. The mask is dropped before the selected cloud's attributes are copied. + +The implementation is +`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-features/src/iss.rs`. +It reuses `spatialrust-search::parallel_worker_count`: one worker below 4,096 +points; otherwise at most available CPU parallelism and +`max(1, floor(points / 16,384))` useful chunks. With default chunking, fewer than +32,768 points retain serial execution. One-worker calls run on the caller and +collect indices directly, without a selection mask. Worker injection is private +and used by tests; no public worker setting, dependency or feature flag is added. + +Unsorted radius-neighbor order, f64 centroid/covariance arithmetic, eigenvalue +sorting, ratio comparisons and f32 saliency conversion retain their original +operation order. There is no floating-point reduction across workers. Non-salient +NaN rows are skipped; equal local maxima all survive the existing `>=` rule. +The old comment describing strict maxima is corrected to match that behavior. +Radius validation, unusual/NaN configurations and empty-cloud behavior remain +unchanged. Selected indices, schema, metadata and all six buffer types are +preserved. The existing covariance/eigenvalue and attribute-gather helpers are +verified unchanged against the parent source. + +Saliency retains four payload bytes per input point. Parallel selection adds +one temporary byte per point; it does not allocate per-worker index vectors +or a second result-sized index merge buffer. KD-tree, selected indices and the +returned cloud remain input/result-sized. Each worker holds only its current +radius-neighbor vector, which may include the entire cloud for a large radius. +Worker limits are per call rather than a process-wide CPU or scratch-byte budget. +Concurrent callers have independent budgets. No all-query neighbor table, +distance matrix, global mutable state or device transfer is introduced. + +## Correctness and feature validation + +The frozen parent module is +`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-features/tests/fixtures/iss_serial.rs`. +It matches source e5b0c6a after line-ending/trailing-newline normalization, with +only the original test module omitted. It is compiled solely in native test +builds and retains the independent serial public detector. + +Four tests in +`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-features/src/iss_parallel_tests.rs` +cover: + +- **1,081,992 exact selection decisions and selected clouds** against the frozen + serial public API: 0/1/17/4,095/4,096/4,103/32,771 points, planar/noisy curved/ + duplicate noisy curved geometry, radii 0.075/0.16, forced 1/2/6 workers and + public staging. Large cases cross real multi-worker staging with uneven bands. +- Independent covariance eigenvalues for eight anisotropic box corners. All + eight equal maxima survive suppression and an isolated ninth point is + excluded. Coincident, collinear and sparse support yield no keypoints. +- Exact selected indices and gathered F32/F64/U8/U16/U32/I32 fields on a + 32,771-point cloud, with renamed XYZ semantics, metadata, signed-zero words + and NaN payload words. Both forced-six-worker and public paths are checked. +- Frozen results/errors for invalid, NaN and infinite radii, unusual/NaN gamma + thresholds, zero/oversized minimum support and empty/singleton input. Zero/ + oversized private worker requests stay safe. Concurrent calls preserve + independent results on shared borrowed input. + +Minimal `feature-iss` tests pass: seven native tests and one ignored benchmark. +All-feature tests pass: twenty-three native tests and two ignored benchmarks, +including boundary detection, CPU/GPU normals and orientation. One facade +smoke passes. Minimal/all-feature all-target strict Clippy and facade smoke +Clippy pass. Crate and facade WASM compilation passes; browser execution or +threading performance is not established. Formatting and diff checks pass. +No workflow or check names change; all 114 named CI checks gate merge. + +```text +cargo test -p spatialrust-features --no-default-features --features feature-iss --lib -- --test-threads=1 --nocapture +cargo test -p spatialrust-features --all-features --tests -- --test-threads=1 +cargo test -p spatialrust --no-default-features --features feature-iss --test iss_parallel_smoke +cargo clippy -p spatialrust-features --no-default-features --features feature-iss --all-targets -- -D warnings +cargo clippy -p spatialrust-features --all-features --all-targets -- -D warnings +cargo clippy -p spatialrust --no-default-features --features feature-iss --test iss_parallel_smoke -- -D warnings +cargo check -p spatialrust-features --target wasm32-unknown-unknown --no-default-features --features feature-iss +cargo check -p spatialrust --target wasm32-unknown-unknown --no-default-features --features feature-iss +cargo fmt --all --check +git diff --check +``` + +## Completed detection-API measurement + +Release native Windows, Rust/Cargo 1.98.1. Twelve conditions combine +512/16,384/32,768/65,536 points with planar, noisy curved and duplicate noisy +curved inputs. Radius is 0.11 for the plane and 0.16 for the other fixtures; +suppression radius is 0.75 times saliency radius. Gamma thresholds and minimum +support retain their defaults. This machine reports 12 available CPU workers; +default policy uses 1/1/2/4 workers respectively. No other tests or builds +launched by this task overlap the run. Desktop activity and other callers' +CPU usage are not controlled. + +The ignored diagnostic is +`C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-features/src/iss_parallel_benchmark.rs`: + +```text +cargo test -p spatialrust-features --no-default-features --features feature-iss --release --lib paired_previous_serial_parallel_iss -- --ignored --nocapture --test-threads=1 +``` + +Each condition runs two warmups and twenty rotating rounds of the frozen serial +public API, the new forced-one-worker path and the current public API. All 792 +completed warmup/measured results are checked against frozen selected indices +and clouds. Timers include validation, KD-tree construction, saliency allocation, +covariance/eigenvalue computation, local-maximum selection, worker spawn/join, +ordered index collection, schema/attribute gathering and scratch/tree destruction. +Fixture construction, output verification and returned result destruction are +outside the timer. These are complete host-returning `detect` measurements, +including the three XYZ attribute columns in these fixtures. + +Table times are per-variant medians. Ratios are medians of twenty paired old/new +times, with 10,000 paired median bootstrap resamples (seed 20261012, sorted +indices 249 and 9749). Above one favors the new path. Intervals describe this +run's samples, rather than portability or repeat-run guarantees. + +| Points | Geometry | Workers | Old / forced one / public ms | Old / forced one (95% CI) | Old / public (95% CI) | +| ---: | --- | ---: | ---: | ---: | ---: | +| 512 | plane | 1 | 0.896 / 0.921 / 0.869 | 1.039 (0.973–1.168) | 0.962 (0.881–1.065) | +| 512 | noisy curved | 1 | 1.497 / 1.560 / 1.459 | 0.987 (0.939–1.094) | 1.006 (0.951–1.058) | +| 512 | duplicate noisy curved | 1 | 1.004 / 0.928 / 0.951 | 1.037 (0.985–1.149) | 0.999 (0.982–1.030) | +| 16,384 | plane | 1 | 36.406 / 35.885 / 36.119 | 0.989 (0.961–1.050) | 1.000 (0.955–1.034) | +| 16,384 | noisy curved | 1 | 52.939 / 53.267 / 59.313 | 0.998 (0.918–1.114) | 0.957 (0.887–0.994) | +| 16,384 | duplicate noisy curved | 1 | 65.307 / 61.678 / 61.220 | 1.031 (0.991–1.129) | 1.013 (0.911–1.070) | +| 32,768 | plane | 2 | 42.016 / 41.405 / 29.342 | 0.953 (0.901–1.088) | 1.411 (1.262–1.577) | +| 32,768 | noisy curved | 2 | 99.390 / 104.419 / 69.727 | 0.902 (0.866–0.956) | 1.450 (1.373–1.492) | +| 32,768 | duplicate noisy curved | 2 | 165.468 / 144.637 / 90.523 | 1.049 (0.949–1.207) | 1.837 (1.431–2.005) | +| 65,536 | plane | 4 | 70.653 / 73.114 / 34.562 | 1.041 (0.918–1.077) | 2.096 (1.946–2.266) | +| 65,536 | noisy curved | 4 | 201.126 / 212.835 / 77.801 | 0.985 (0.865–1.074) | 2.541 (2.382–2.729) | +| 65,536 | duplicate noisy curved | 4 | 299.961 / 312.491 / 117.359 | 0.918 (0.827–1.023) | 2.321 (2.113–2.670) | + +All six two/four-worker conditions favor public staging, with paired median +speedups 1.411–2.541 and interval lower bounds above one. The 16,384-point noisy +curved public control has an adverse ratio of 0.957 (0.887–0.994); the other +five one-worker public intervals include one. The forced-one-worker noisy +curved case at 32,768 points also regresses, at 0.902 (0.866–0.956). Both adverse +controls are retained; serial neutrality and fixed cross-machine speedups are +not claimed. Timing variation is visible throughout this uncontrolled desktop +run, and bootstrap intervals do not capture all sources of systematic variation. +The existing conservative worker policy is retained. All measured public +parallel paths complete the same detection API, including selected XYZ copying; +there is no claim about GPU ISS or a whole registration pipeline. + +Raw samples and analysis remain private ignored files under +`C:/Users/rsasa/Workspace/SpatialRust/target/automation/parallel-iss-*`. +Published measurements describe execution conditions without a local adapter +model name. GPU ISS, parallel attribute gathering and process-wide worker +budgets are outside this slice.