spatialrust_features/
boundary.rs

1//! Boundary / edge point detection via tangent-plane angle gaps.
2//!
3//! A point on the interior of a well-sampled surface has neighbors spread all
4//! around it in its tangent plane; a point on a boundary (a hole rim or the edge
5//! of a scan) has a large empty angular wedge. Projecting neighbors into the
6//! tangent plane and measuring the largest gap between consecutive directions
7//! flags those boundary points.
8
9use spatialrust_core::{
10    HasNormals3, HasPositions3, PointBuffer, PointBufferSet, PointCloud, SpatialError,
11    SpatialResult,
12};
13use spatialrust_math::Vec3;
14use spatialrust_search::{KdTree, RadiusSearchIndex};
15
16/// Configuration for [`BoundaryDetector`].
17#[derive(Clone, Copy, Debug, PartialEq)]
18pub struct BoundaryConfig {
19    /// Neighborhood radius used to gather tangent-plane directions.
20    pub search_radius: f32,
21    /// A point is a boundary if its largest angular gap (radians) exceeds this.
22    pub angle_threshold: f32,
23    /// Minimum neighbors required to evaluate a point.
24    pub min_neighbors: usize,
25}
26
27impl Default for BoundaryConfig {
28    fn default() -> Self {
29        // ~90°: interior points have small gaps, boundary points a large wedge.
30        Self { search_radius: 0.1, angle_threshold: std::f32::consts::FRAC_PI_2, min_neighbors: 5 }
31    }
32}
33
34impl BoundaryConfig {
35    /// Creates a config from the search radius (90° threshold).
36    #[must_use]
37    pub fn with_radius(search_radius: f32) -> Self {
38        Self { search_radius, ..Self::default() }
39    }
40}
41
42/// Result of boundary detection.
43#[derive(Clone, Debug, PartialEq)]
44pub struct BoundaryResult {
45    /// Indices of the boundary points in the input cloud.
46    pub indices: Vec<usize>,
47    /// Sub-cloud containing only the boundary points (schema preserved).
48    pub boundary: PointCloud,
49}
50
51/// Tangent-plane boundary point detector. The input cloud must carry normals.
52#[derive(Clone, Copy, Debug, PartialEq)]
53pub struct BoundaryDetector {
54    config: BoundaryConfig,
55}
56
57impl BoundaryDetector {
58    /// Creates a detector from config.
59    #[must_use]
60    pub const fn new(config: BoundaryConfig) -> Self {
61        Self { config }
62    }
63
64    /// Returns the detector config.
65    #[must_use]
66    pub const fn config(&self) -> BoundaryConfig {
67        self.config
68    }
69
70    /// Computes the per-point boundary mask (`true` = boundary).
71    pub fn boundary_mask(&self, input: &PointCloud) -> SpatialResult<Vec<bool>> {
72        if self.config.search_radius <= 0.0 || self.config.search_radius.is_nan() {
73            return Err(SpatialError::InvalidArgument("search_radius must be positive".to_owned()));
74        }
75        let len = input.len();
76        if len == 0 {
77            return Ok(Vec::new());
78        }
79
80        let (x, y, z) = input.positions3()?;
81        let (nx, ny, nz) = input.normals3()?;
82        let tree = KdTree::from_slices(x, y, z);
83
84        let mut mask = vec![false; len];
85        let mut angles: Vec<f32> = Vec::new();
86        for i in 0..len {
87            let p = Vec3::new(x[i], y[i], z[i]);
88            let normal = Vec3::new(nx[i], ny[i], nz[i]);
89            let (u, v) = tangent_basis(normal);
90
91            let neighbors = tree.radius_search(p.x, p.y, p.z, self.config.search_radius);
92            angles.clear();
93            for neighbor in &neighbors {
94                let j = neighbor.index;
95                if j == i {
96                    continue;
97                }
98                let d = Vec3::new(x[j] - p.x, y[j] - p.y, z[j] - p.z);
99                let (du, dv) = (d.dot(u), d.dot(v));
100                if du * du + dv * dv > 1e-12 {
101                    angles.push(dv.atan2(du));
102                }
103            }
104            if angles.len() < self.config.min_neighbors {
105                continue;
106            }
107            mask[i] = max_angle_gap(&mut angles) > self.config.angle_threshold;
108        }
109        Ok(mask)
110    }
111
112    /// Detects boundary points, returning their indices and a sub-cloud.
113    pub fn detect(&self, input: &PointCloud) -> SpatialResult<BoundaryResult> {
114        let mask = self.boundary_mask(input)?;
115        let indices: Vec<usize> =
116            mask.iter().enumerate().filter_map(|(i, &b)| b.then_some(i)).collect();
117        let boundary = gather_indices(input, &indices)?;
118        Ok(BoundaryResult { indices, boundary })
119    }
120}
121
122/// Builds an orthonormal basis `(u, v)` spanning the plane perpendicular to `n`.
123fn tangent_basis(n: Vec3<f32>) -> (Vec3<f32>, Vec3<f32>) {
124    let helper = if n.x.abs() < 0.9 { Vec3::new(1.0, 0.0, 0.0) } else { Vec3::new(0.0, 1.0, 0.0) };
125    let u = n.cross(helper).normalize();
126    let v = n.cross(u);
127    (u, v)
128}
129
130/// Largest gap between consecutive sorted angles, wrapping around the circle.
131fn max_angle_gap(angles: &mut [f32]) -> f32 {
132    // A point with no in-plane neighbors has no gap to measure (reachable when
133    // `min_neighbors` is 0 and a point is isolated); treat it as non-salient.
134    if angles.is_empty() {
135        return 0.0;
136    }
137    angles.sort_by(f32::total_cmp);
138    let mut max_gap = 0.0_f32;
139    for w in angles.windows(2) {
140        max_gap = max_gap.max(w[1] - w[0]);
141    }
142    // Wrap-around gap between the last and first direction.
143    let wrap = (angles[0] + std::f32::consts::TAU) - angles[angles.len() - 1];
144    max_gap.max(wrap)
145}
146
147fn gather_indices(input: &PointCloud, indices: &[usize]) -> SpatialResult<PointCloud> {
148    let mut buffers = PointBufferSet::new();
149    for field in input.schema().fields() {
150        let source = input.field(&field.name)?;
151        buffers.insert(field.name.clone(), gather_buffer(source, indices));
152    }
153    PointCloud::try_from_parts(input.schema().clone(), buffers, input.metadata().clone())
154}
155
156fn gather_buffer(source: &PointBuffer, indices: &[usize]) -> PointBuffer {
157    match source {
158        PointBuffer::F32(v) => PointBuffer::from_f32(indices.iter().map(|&i| v[i]).collect()),
159        PointBuffer::F64(v) => PointBuffer::F64(indices.iter().map(|&i| v[i]).collect()),
160        PointBuffer::U8(v) => PointBuffer::U8(indices.iter().map(|&i| v[i]).collect()),
161        PointBuffer::U16(v) => PointBuffer::U16(indices.iter().map(|&i| v[i]).collect()),
162        PointBuffer::U32(v) => PointBuffer::U32(indices.iter().map(|&i| v[i]).collect()),
163        PointBuffer::I32(v) => PointBuffer::I32(indices.iter().map(|&i| v[i]).collect()),
164    }
165}
166
167#[cfg(test)]
168mod tests {
169    use super::{BoundaryConfig, BoundaryDetector};
170    use spatialrust_core::{DType, FieldSemantic, PointCloudBuilder, PointField, PointSchema};
171
172    fn schema() -> PointSchema {
173        PointSchema::new()
174            .with_field(PointField::scalar("x", FieldSemantic::PositionX, DType::F32))
175            .with_field(PointField::scalar("y", FieldSemantic::PositionY, DType::F32))
176            .with_field(PointField::scalar("z", FieldSemantic::PositionZ, DType::F32))
177            .with_field(PointField::scalar("normal_x", FieldSemantic::NormalX, DType::F32))
178            .with_field(PointField::scalar("normal_y", FieldSemantic::NormalY, DType::F32))
179            .with_field(PointField::scalar("normal_z", FieldSemantic::NormalZ, DType::F32))
180    }
181
182    /// A square planar patch (normal +Z): rim points are boundaries, the
183    /// interior is not.
184    fn square_patch(n: usize) -> spatialrust_core::PointCloud {
185        let mut builder = PointCloudBuilder::new(schema());
186        for i in 0..n {
187            for j in 0..n {
188                builder.push_point([i as f32 * 0.1, j as f32 * 0.1, 0.0, 0.0, 0.0, 1.0]).unwrap();
189            }
190        }
191        builder.build().unwrap()
192    }
193
194    #[test]
195    fn detects_patch_rim_not_interior() {
196        let n = 12;
197        let cloud = square_patch(n);
198        // min_neighbors=3 so sparse corners (only 3 in-radius neighbors) are
199        // still evaluated.
200        let config = BoundaryConfig { min_neighbors: 3, ..BoundaryConfig::with_radius(0.18) };
201        let result = BoundaryDetector::new(config).detect(&cloud).unwrap();
202        let mask = BoundaryDetector::new(config).boundary_mask(&cloud).unwrap();
203
204        // The center point is interior; a corner is on the boundary.
205        let center = (n / 2) * n + (n / 2);
206        let corner = 0;
207        assert!(!mask[center], "interior point flagged as boundary");
208        assert!(mask[corner], "corner not flagged as boundary");
209        // Boundaries are a minority of the patch.
210        assert!(!result.indices.is_empty());
211        assert!(result.indices.len() < cloud.len() / 2);
212        assert_eq!(result.boundary.len(), result.indices.len());
213    }
214
215    #[test]
216    fn rejects_bad_radius() {
217        let cloud = square_patch(5);
218        assert!(BoundaryDetector::new(BoundaryConfig::with_radius(0.0)).detect(&cloud).is_err());
219    }
220
221    #[test]
222    fn handles_isolated_point_with_zero_min_neighbors() {
223        // A lone point far from the patch has no in-radius neighbors. With
224        // min_neighbors = 0 this used to reach the angle-gap code with an empty
225        // angle list and panic; it must now just run.
226        let mut builder = PointCloudBuilder::new(schema());
227        builder.push_point([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]).unwrap();
228        builder.push_point([100.0, 100.0, 100.0, 0.0, 0.0, 1.0]).unwrap();
229        let cloud = builder.build().unwrap();
230
231        let config = BoundaryConfig { min_neighbors: 0, ..BoundaryConfig::with_radius(0.1) };
232        let mask = BoundaryDetector::new(config).boundary_mask(&cloud).unwrap();
233        // Isolated points (no neighbors) are not flagged as surface boundaries.
234        assert_eq!(mask, vec![false, false]);
235    }
236}
spatialrust_features/boundary.rs

spatialrust_features/
boundary.rs
