brepkit_geometry/sampling/
surface.rs

1//! Regular grid sampling of parametric surfaces.
2
3use brepkit_math::traits::ParametricSurface;
4use brepkit_math::vec::Point3;
5
6/// Sample a regular N×M grid of surface points.
7///
8/// Returns `nu` rows, each with `nv` points.
9///
10/// - Row `i` corresponds to `u = u_range.0 + i*(u_range.1 - u_range.0) / (nu - 1)`.
11/// - Column `j` corresponds to `v = v_range.0 + j*(v_range.1 - v_range.0) / (nv - 1)`.
12///
13/// Edge cases:
14/// - If `nu == 0` or `nv == 0`, an empty `Vec` is returned.
15/// - If `nu == 1`, the single row is evaluated at `u_range.0`.
16/// - If `nv == 1`, each row contains one point evaluated at `v_range.0`.
17#[must_use]
18pub fn surface_grid<S: ParametricSurface>(
19    surface: &S,
20    u_range: (f64, f64),
21    v_range: (f64, f64),
22    nu: usize,
23    nv: usize,
24) -> Vec<Vec<Point3>> {
25    if nu == 0 || nv == 0 {
26        return Vec::new();
27    }
28
29    (0..nu)
30        .map(|i| {
31            let u = if nu == 1 {
32                u_range.0
33            } else if i == nu - 1 {
34                u_range.1
35            } else {
36                u_range.0 + i as f64 * (u_range.1 - u_range.0) / (nu - 1) as f64
37            };
38
39            (0..nv)
40                .map(|j| {
41                    let v = if nv == 1 {
42                        v_range.0
43                    } else if j == nv - 1 {
44                        v_range.1
45                    } else {
46                        v_range.0 + j as f64 * (v_range.1 - v_range.0) / (nv - 1) as f64
47                    };
48                    surface.evaluate(u, v)
49                })
50                .collect()
51        })
52        .collect()
53}
54
55#[cfg(test)]
56mod tests {
57    #![allow(clippy::unwrap_used, clippy::expect_used)]
58
59    use super::*;
60    use brepkit_math::surfaces::CylindricalSurface;
61    use brepkit_math::vec::{Point3, Vec3};
62    use std::f64::consts::TAU;
63
64    fn unit_cylinder() -> CylindricalSurface {
65        CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0).unwrap()
66    }
67
68    #[test]
69    fn zero_nu_returns_empty() {
70        let s = unit_cylinder();
71        let grid = surface_grid(&s, (0.0, TAU), (0.0, 1.0), 0, 4);
72        assert!(grid.is_empty());
73    }
74
75    #[test]
76    fn zero_nv_returns_empty() {
77        let s = unit_cylinder();
78        let grid = surface_grid(&s, (0.0, TAU), (0.0, 1.0), 4, 0);
79        assert!(grid.is_empty());
80    }
81
82    #[test]
83    fn grid_dimensions_correct() {
84        let s = unit_cylinder();
85        let grid = surface_grid(&s, (0.0, TAU), (0.0, 1.0), 5, 3);
86        assert_eq!(grid.len(), 5, "expected 5 rows");
87        for (i, row) in grid.iter().enumerate() {
88            assert_eq!(row.len(), 3, "row {i} should have 3 columns");
89        }
90    }
91
92    #[test]
93    fn all_points_on_cylinder_surface() {
94        // Cylinder P(u,v) = (cos u, sin u, v).  All points must satisfy x²+y²=1.
95        let s = unit_cylinder();
96        let grid = surface_grid(&s, (0.0, TAU), (0.0, 2.0), 8, 4);
97        for row in &grid {
98            for p in row {
99                let r = (p.x() * p.x() + p.y() * p.y()).sqrt();
100                assert!((r - 1.0).abs() < 1e-10, "point not on cylinder: r={r}");
101            }
102        }
103    }
104
105    #[test]
106    fn first_and_last_u_are_endpoints() {
107        let s = unit_cylinder();
108        let u0 = 0.0_f64;
109        let u1 = TAU;
110        let grid = surface_grid(&s, (u0, u1), (0.0, 1.0), 4, 2);
111        assert_eq!(grid.len(), 4);
112
113        // First row: u = u0 → point should equal surface.evaluate(u0, 0.0)
114        let expected_first = s.evaluate(u0, 0.0);
115        let actual_first = grid[0][0];
116        let d0 = {
117            let dx = actual_first.x() - expected_first.x();
118            let dy = actual_first.y() - expected_first.y();
119            let dz = actual_first.z() - expected_first.z();
120            (dx * dx + dy * dy + dz * dz).sqrt()
121        };
122        assert!(d0 < 1e-12, "first row u mismatch: dist={d0}");
123
124        // Last row: u = u1.
125        let expected_last = s.evaluate(u1, 0.0);
126        let actual_last = grid[3][0];
127        let d1 = {
128            let dx = actual_last.x() - expected_last.x();
129            let dy = actual_last.y() - expected_last.y();
130            let dz = actual_last.z() - expected_last.z();
131            (dx * dx + dy * dy + dz * dz).sqrt()
132        };
133        assert!(d1 < 1e-12, "last row u mismatch: dist={d1}");
134    }
135
136    #[test]
137    fn single_row_and_column() {
138        let s = unit_cylinder();
139        let grid = surface_grid(&s, (0.0, TAU), (0.5, 1.5), 1, 1);
140        assert_eq!(grid.len(), 1);
141        assert_eq!(grid[0].len(), 1);
142        let expected = s.evaluate(0.0, 0.5);
143        let p = grid[0][0];
144        let d = {
145            let dx = p.x() - expected.x();
146            let dy = p.y() - expected.y();
147            let dz = p.z() - expected.z();
148            (dx * dx + dy * dy + dz * dz).sqrt()
149        };
150        assert!(d < 1e-12);
151    }
152
153    #[test]
154    fn v_endpoint_included_in_last_column() {
155        let s = unit_cylinder();
156        let v0 = 0.0_f64;
157        let v1 = 3.0_f64;
158        let grid = surface_grid(&s, (0.0, 1.0), (v0, v1), 2, 4);
159        // Last column (j=3) should have v = v1.
160        for (i, row) in grid.iter().enumerate() {
161            let expected = s.evaluate(if i == 0 { 0.0 } else { 1.0 }, v1);
162            let actual = row[3];
163            let d = {
164                let dx = actual.x() - expected.x();
165                let dy = actual.y() - expected.y();
166                let dz = actual.z() - expected.z();
167                (dx * dx + dy * dy + dz * dz).sqrt()
168            };
169            assert!(d < 1e-12, "v endpoint mismatch at row {i}: dist={d}");
170        }
171    }
172}
brepkit_geometry/sampling/surface.rs

brepkit_geometry/sampling/
surface.rs
