spatialrust_features/
iss.rs

1//! Intrinsic Shape Signatures (ISS) keypoint detection.
2//!
3//! ISS picks a sparse set of geometrically salient points — corners and other
4//! spots with variation in all three directions — by thresholding the ratios of
5//! the eigenvalues of each point's neighborhood covariance, then keeping local
6//! maxima of the smallest eigenvalue. The keypoints are a natural front-end for
7//! descriptor-based global registration: far fewer points to describe and match.
8
9use spatialrust_core::{
10    HasPositions3, PointBuffer, PointBufferSet, PointCloud, SpatialError, SpatialResult,
11};
12use spatialrust_math::{symmetric_eigen3, Mat3};
13use spatialrust_search::{KdTree, RadiusSearchIndex};
14
15/// Configuration for [`IssKeypointDetector`].
16#[derive(Clone, Copy, Debug, PartialEq)]
17pub struct IssKeypointConfig {
18    /// Radius of the neighborhood used to build each covariance matrix.
19    pub salient_radius: f32,
20    /// Radius for non-maximum suppression of the saliency (smallest eigenvalue).
21    pub non_max_radius: f32,
22    /// Upper bound on `lambda2 / lambda1` (rejects flat / planar points).
23    pub gamma_21: f32,
24    /// Upper bound on `lambda3 / lambda2` (rejects edge / linear points).
25    pub gamma_32: f32,
26    /// Minimum neighbors within `salient_radius` for a point to be considered.
27    pub min_neighbors: usize,
28}
29
30impl Default for IssKeypointConfig {
31    fn default() -> Self {
32        Self {
33            salient_radius: 0.2,
34            non_max_radius: 0.15,
35            gamma_21: 0.975,
36            gamma_32: 0.975,
37            min_neighbors: 5,
38        }
39    }
40}
41
42impl IssKeypointConfig {
43    /// Creates a config from the salient and non-max-suppression radii.
44    #[must_use]
45    pub fn with_radii(salient_radius: f32, non_max_radius: f32) -> Self {
46        Self { salient_radius, non_max_radius, ..Self::default() }
47    }
48}
49
50/// Result of ISS keypoint detection.
51#[derive(Clone, Debug, PartialEq)]
52pub struct IssKeypointResult {
53    /// Indices of the keypoints in the input cloud.
54    pub indices: Vec<usize>,
55    /// Sub-cloud containing only the keypoints (input schema preserved).
56    pub keypoints: PointCloud,
57}
58
59/// Intrinsic Shape Signatures keypoint detector.
60#[derive(Clone, Copy, Debug, PartialEq)]
61pub struct IssKeypointDetector {
62    config: IssKeypointConfig,
63}
64
65impl IssKeypointDetector {
66    /// Creates a detector from config.
67    #[must_use]
68    pub const fn new(config: IssKeypointConfig) -> Self {
69        Self { config }
70    }
71
72    /// Returns the detector config.
73    #[must_use]
74    pub const fn config(&self) -> IssKeypointConfig {
75        self.config
76    }
77
78    /// Detects keypoints and returns their indices and a keypoint sub-cloud.
79    pub fn detect(&self, input: &PointCloud) -> SpatialResult<IssKeypointResult> {
80        if self.config.salient_radius <= 0.0 || self.config.non_max_radius <= 0.0 {
81            return Err(SpatialError::InvalidArgument("ISS radii must be positive".to_owned()));
82        }
83
84        let (x, y, z) = input.positions3()?;
85        let len = input.len();
86        let tree = KdTree::from_slices(x, y, z);
87
88        // Saliency = smallest eigenvalue (lambda3) where the eigenvalue-ratio
89        // tests pass; NaN marks a non-salient point.
90        let mut saliency = vec![f32::NAN; len];
91        for i in 0..len {
92            let neighbors = tree.radius_search(x[i], y[i], z[i], self.config.salient_radius);
93            if neighbors.len() < self.config.min_neighbors {
94                continue;
95            }
96            let Some(eigenvalues) = neighborhood_eigenvalues(x, y, z, &neighbors) else {
97                continue;
98            };
99            // Ascending order: l3 <= l2 <= l1.
100            let (l3, l2, l1) = (eigenvalues[0], eigenvalues[1], eigenvalues[2]);
101            if l1 <= 0.0 || l2 <= 0.0 {
102                continue;
103            }
104            if (l2 / l1) < f64::from(self.config.gamma_21)
105                && (l3 / l2) < f64::from(self.config.gamma_32)
106            {
107                saliency[i] = l3 as f32;
108            }
109        }
110
111        // Non-maximum suppression: keep points whose saliency is the strict
112        // maximum within `non_max_radius`.
113        let mut indices = Vec::new();
114        for i in 0..len {
115            let s = saliency[i];
116            if s.is_nan() {
117                continue;
118            }
119            let neighbors = tree.radius_search(x[i], y[i], z[i], self.config.non_max_radius);
120            let is_local_max = neighbors.iter().all(|n| {
121                let other = saliency[n.index];
122                n.index == i || other.is_nan() || s >= other
123            });
124            if is_local_max {
125                indices.push(i);
126            }
127        }
128
129        let keypoints = gather_indices(input, &indices)?;
130        Ok(IssKeypointResult { indices, keypoints })
131    }
132}
133
134/// Eigenvalues (ascending) of the neighborhood covariance about its centroid.
135fn neighborhood_eigenvalues(
136    x: &[f32],
137    y: &[f32],
138    z: &[f32],
139    neighbors: &[spatialrust_search::Neighbor],
140) -> Option<[f64; 3]> {
141    let count = neighbors.len() as f64;
142    if count < 3.0 {
143        return None;
144    }
145    let (mut mx, mut my, mut mz) = (0.0_f64, 0.0_f64, 0.0_f64);
146    for n in neighbors {
147        mx += f64::from(x[n.index]);
148        my += f64::from(y[n.index]);
149        mz += f64::from(z[n.index]);
150    }
151    mx /= count;
152    my /= count;
153    mz /= count;
154
155    let (mut c00, mut c11, mut c22) = (0.0_f64, 0.0, 0.0);
156    let (mut c01, mut c02, mut c12) = (0.0_f64, 0.0, 0.0);
157    for n in neighbors {
158        let dx = f64::from(x[n.index]) - mx;
159        let dy = f64::from(y[n.index]) - my;
160        let dz = f64::from(z[n.index]) - mz;
161        c00 += dx * dx;
162        c11 += dy * dy;
163        c22 += dz * dz;
164        c01 += dx * dy;
165        c02 += dx * dz;
166        c12 += dy * dz;
167    }
168    let inv = 1.0 / count;
169    let covariance = Mat3::<f64>::from_rows(
170        [c00 * inv, c01 * inv, c02 * inv],
171        [c01 * inv, c11 * inv, c12 * inv],
172        [c02 * inv, c12 * inv, c22 * inv],
173    );
174    Some(symmetric_eigen3(covariance).eigenvalues)
175}
176
177/// Gathers the selected indices into a new cloud, preserving schema.
178fn gather_indices(input: &PointCloud, indices: &[usize]) -> SpatialResult<PointCloud> {
179    let mut buffers = PointBufferSet::new();
180    for field in input.schema().fields() {
181        let source = input.field(&field.name)?;
182        buffers.insert(field.name.clone(), gather_buffer(source, indices));
183    }
184    PointCloud::try_from_parts(input.schema().clone(), buffers, input.metadata().clone())
185}
186
187fn gather_buffer(source: &PointBuffer, indices: &[usize]) -> PointBuffer {
188    match source {
189        PointBuffer::F32(v) => PointBuffer::from_f32(indices.iter().map(|&i| v[i]).collect()),
190        PointBuffer::F64(v) => PointBuffer::F64(indices.iter().map(|&i| v[i]).collect()),
191        PointBuffer::U8(v) => PointBuffer::U8(indices.iter().map(|&i| v[i]).collect()),
192        PointBuffer::U16(v) => PointBuffer::U16(indices.iter().map(|&i| v[i]).collect()),
193        PointBuffer::U32(v) => PointBuffer::U32(indices.iter().map(|&i| v[i]).collect()),
194        PointBuffer::I32(v) => PointBuffer::I32(indices.iter().map(|&i| v[i]).collect()),
195    }
196}
197
198#[cfg(test)]
199mod tests {
200    use super::{IssKeypointConfig, IssKeypointDetector};
201    use spatialrust_core::{PointCloudBuilder, StandardSchemas};
202
203    /// A flat floor plus a sharp spike — the spike tip should be a keypoint and
204    /// the flat region should not.
205    fn floor_with_spike() -> spatialrust_core::PointCloud {
206        let mut builder = PointCloudBuilder::new(StandardSchemas::point_xyz());
207        for i in 0..30 {
208            for j in 0..30 {
209                builder.push_point([i as f32 * 0.05, j as f32 * 0.05, 0.0]).unwrap();
210            }
211        }
212        // A small dense cluster rising off the plane near its center.
213        for k in 0..40 {
214            let h = k as f32 * 0.02;
215            builder.push_point([0.75, 0.75, h]).unwrap();
216            builder.push_point([0.77, 0.75, h]).unwrap();
217            builder.push_point([0.75, 0.77, h]).unwrap();
218        }
219        builder.build().unwrap()
220    }
221
222    #[test]
223    fn finds_keypoints_on_salient_structure() {
224        let cloud = floor_with_spike();
225        let detector = IssKeypointDetector::new(IssKeypointConfig {
226            salient_radius: 0.12,
227            non_max_radius: 0.08,
228            gamma_21: 0.9,
229            gamma_32: 0.9,
230            min_neighbors: 5,
231        });
232        let result = detector.detect(&cloud).unwrap();
233        // Some keypoints found, far fewer than the input size.
234        assert!(!result.indices.is_empty());
235        assert!(result.indices.len() < cloud.len() / 4);
236        assert_eq!(result.keypoints.len(), result.indices.len());
237    }
238
239    #[test]
240    fn flat_plane_interior_is_not_salient() {
241        let n = 20;
242        let mut builder = PointCloudBuilder::new(StandardSchemas::point_xyz());
243        for i in 0..n {
244            for j in 0..n {
245                builder.push_point([i as f32 * 0.05, j as f32 * 0.05, 0.0]).unwrap();
246            }
247        }
248        let cloud = builder.build().unwrap();
249        let detector = IssKeypointDetector::new(IssKeypointConfig::with_radii(0.12, 0.08));
250        let result = detector.detect(&cloud).unwrap();
251        // ISS legitimately flags the plane's boundary (a real edge), but the
252        // isotropic interior (lambda1 ~ lambda2) must be rejected by gamma_21.
253        let center = (n / 2) * n + (n / 2);
254        assert!(!result.indices.contains(&center), "plane interior should not be salient");
255    }
256
257    #[test]
258    fn rejects_bad_radii() {
259        let cloud = floor_with_spike();
260        let detector = IssKeypointDetector::new(IssKeypointConfig::with_radii(0.0, 0.1));
261        assert!(detector.detect(&cloud).is_err());
262    }
263}
spatialrust_features/iss.rs

spatialrust_features/
iss.rs
