geop_core_geometry/nurb_surface/
curvature.rs

1use geop_core_math::{geop_error::GeopResult, scalars::Scalar};
2
3use super::NurbSurface;
4
5impl<S: Scalar> NurbSurface<S, 4> {
6    /// A conservative radius of curvature at `(u, v)`, or `None` if the
7    /// surface is (locally) flat in both parametric directions.
8    ///
9    /// Computed straight from the surface's own second partial derivatives
10    /// — no history, no finite differences — via the normal curvature along
11    /// each parametric direction, `κ_a = (S_aa · n) / |S_a|²` for `a ∈ {u,
12    /// v}` (the diagonal terms of the second fundamental form divided by the
13    /// diagonal terms of the first). This is the *exact* normal curvature in
14    /// direction `a` only when `Su ⊥ Sv`; the general formula also needs the
15    /// mixed partial `Suv` and the off-diagonal metric term `F = Su·Sv` to
16    /// handle an arbitrary direction. Every surface this crate actually
17    /// constructs has orthogonal parametric directions by construction —
18    /// flat bilinear box/cap faces (`Su`, `Sv` are the patch's two edge
19    /// directions) and `revolve`'s ruled patches (axial `u` is always
20    /// perpendicular to the circular `v`) — so the approximation is exact
21    /// for our surfaces, not just a rough heuristic.
22    ///
23    /// The returned radius is `1 / max(|κ_u|, |κ_v|)`: the tighter of the
24    /// two bends, so a caller sizing steps off of it stays conservative.
25    pub fn curvature_radius(&self, u: S, v: S) -> GeopResult<Option<S>> {
26        let (su, sv) = self.derivatives(u, v)?;
27        let (suu, svv) = self.second_derivatives(u, v)?;
28
29        let normal = match su.prod_cross(&sv).normalize() {
30            Ok(n) => n,
31            Err(_) => return Ok(None),
32        };
33
34        let su_len_sq = su.prod_dot(&su);
35        let sv_len_sq = sv.prod_dot(&sv);
36
37        let kappa_u = if su_len_sq.could_be_equal(S::ZERO) {
38            S::ZERO
39        } else {
40            suu.prod_dot(&normal).div(su_len_sq)?.abs()
41        };
42        let kappa_v = if sv_len_sq.could_be_equal(S::ZERO) {
43            S::ZERO
44        } else {
45            svv.prod_dot(&normal).div(sv_len_sq)?.abs()
46        };
47
48        let kappa = if kappa_u.definitely_greater(kappa_v) {
49            kappa_u
50        } else {
51            kappa_v
52        };
53
54        if kappa.could_be_equal(S::ZERO) {
55            Ok(None)
56        } else {
57            Ok(Some(S::ONE.div(kappa)?.abs()))
58        }
59    }
60}
61
62#[cfg(test)]
63mod tests {
64    use crate::nurb_surface::NurbSurface;
65    use geop_core_math::for_all_scalars;
66    use geop_core_math::{scalars::Scalar, vector::Vector4};
67
68    fn pt<S: Scalar>(x: f64, y: f64, z: f64, w: f64) -> Vector4<S> {
69        Vector4::from_array([
70            S::from_f64(x),
71            S::from_f64(y),
72            S::from_f64(z),
73            S::from_f64(w),
74        ])
75    }
76
77    /// A flat bilinear patch has zero curvature everywhere: no radius limit.
78    fn check_flat_patch_has_no_curvature_radius<S: Scalar>() {
79        let f = S::from_f64;
80        let s = NurbSurface::try_new(
81            1,
82            1,
83            vec![
84                pt(0., 0., 0., 1.),
85                pt(0., 1., 0., 1.),
86                pt(1., 0., 0., 1.),
87                pt(1., 1., 0., 1.),
88            ],
89            vec![f(0.), f(0.), f(1.), f(1.)],
90            vec![f(0.), f(0.), f(1.), f(1.)],
91        )
92        .unwrap();
93
94        assert!(s.curvature_radius(f(0.5), f(0.5)).unwrap().is_none());
95    }
96    #[test]
97    fn flat_patch_has_no_curvature_radius() {
98        for_all_scalars!(check_flat_patch_has_no_curvature_radius);
99    }
100
101    // A cylindrical-patch regression test (checking `curvature_radius`
102    // against a real `cylinder_solid`) used to live here, commented out.
103    // Since the crate split it can no longer live in this crate at all —
104    // `cylinder_solid` is built by `geop-ops-extrude-revolve`, which sits
105    // *above* `geop-core-geometry` in the workspace's dependency order — so
106    // it was removed rather than carried forward commented out. Equivalent
107    // coverage belongs in `geop-ops-extrude-revolve::cylinder`'s own tests.
108}
geop_core_geometry/nurb_surface/curvature.rs

geop_core_geometry/nurb_surface/
curvature.rs
