Compute GPU ISS saliency and keypoints from resident radius neighborhoods

.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

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

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:

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"]

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<GpuIssKeypointResult> {
+        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,

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};
 

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<u32>,
+    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 }
+}