geop_core_geometry/nurb_surface/
curvature.rs

1use geop_core_math::{
2    geop_error::{GeopError, GeopResult},
3    scalars::Scalar,
4    vector::Vector3,
5};
6
7use super::NurbSurface;
8
9impl NurbSurface<4> {
10    /// A conservative radius of curvature at `(u, v)`, or `None` if the
11    /// surface is (locally) flat in both parametric directions.
12    ///
13    /// Computed straight from the surface's own second partial derivatives
14    /// — no history, no finite differences — via the normal curvature along
15    /// each parametric direction, `κ_a = (S_aa · n) / |S_a|²` for `a ∈ {u,
16    /// v}` (the diagonal terms of the second fundamental form divided by the
17    /// diagonal terms of the first). This is the *exact* normal curvature in
18    /// direction `a` only when `Su ⊥ Sv`; the general formula also needs the
19    /// mixed partial `Suv` and the off-diagonal metric term `F = Su·Sv` to
20    /// handle an arbitrary direction (see [`Self::mean_curvature`]). Every surface this crate actually
21    /// constructs has orthogonal parametric directions by construction —
22    /// flat bilinear box/cap faces (`Su`, `Sv` are the patch's two edge
23    /// directions) and `revolve`'s ruled patches (axial `u` is always
24    /// perpendicular to the circular `v`) — so the approximation is exact
25    /// for our surfaces, not just a rough heuristic.
26    ///
27    /// The returned radius is `1 / max(|κ_u|, |κ_v|)`: the tighter of the
28    /// two bends, so a caller sizing steps off of it stays conservative.
29    pub fn curvature_radius(&self, u: Scalar, v: Scalar) -> GeopResult<Option<Scalar>> {
30        let (su, sv) = self.derivatives(u, v)?;
31        let [suu, _, svv] = self.second_derivatives(u, v)?;
32
33        let normal = match su.prod_cross(&sv).normalize() {
34            Ok(n) => n,
35            Err(_) => return Ok(None),
36        };
37
38        let su_len_sq = su.prod_dot(&su);
39        let sv_len_sq = sv.prod_dot(&sv);
40
41        let kappa_u = if su_len_sq.could_be_equal(Scalar::ZERO) {
42            Scalar::ZERO
43        } else {
44            suu.prod_dot(&normal).div(su_len_sq)?.abs()
45        };
46        let kappa_v = if sv_len_sq.could_be_equal(Scalar::ZERO) {
47            Scalar::ZERO
48        } else {
49            svv.prod_dot(&normal).div(sv_len_sq)?.abs()
50        };
51
52        let kappa = if kappa_u.definitely_greater(kappa_v) {
53            kappa_u
54        } else {
55            kappa_v
56        };
57
58        if kappa.could_be_equal(Scalar::ZERO) {
59            Ok(None)
60        } else {
61            Ok(Some(Scalar::ONE.div(kappa)?.abs()))
62        }
63    }
64
65    /// A radius of curvature at `(u, v)` along the tangent direction
66    /// `along` — a unit vector in the tangent plane, the way a curve on the
67    /// surface runs — or `None` if the surface runs straight that way: the
68    /// normal curvature `κ = (L a² + 2 M a b + N b²)`, where `along = a Su +
69    /// b Sv` and `L, M, N` are the second fundamental form (`Suu·n`,
70    /// `Suv·n`, `Svv·n`). Exact for any parametrization.
71    ///
72    /// A curve on the surface turns, as far as the surface makes it, by this
73    /// — where [`Self::curvature_radius`], the tightest bend in *any*
74    /// direction, says how a curve circling there would. At a cone's apex the
75    /// tightest bend has the radius of the circle round the axis, which goes
76    /// to nothing; a generator running into the apex does not bend at all.
77    ///
78    /// An error where `Su × Sv` vanishes (a pole), and where `along` leaves
79    /// the tangent plane's reach.
80    pub fn curvature_radius_along(
81        &self,
82        u: Scalar,
83        v: Scalar,
84        along: &Vector3,
85    ) -> GeopResult<Option<Scalar>> {
86        let kappa = self.normal_curvature(u, v, along)?.abs();
87        if kappa.could_be_equal(Scalar::ZERO) {
88            Ok(None)
89        } else {
90            Ok(Some(Scalar::ONE.div(kappa)?.abs()))
91        }
92    }
93
94    /// The normal curvature at `(u, v)` along the unit tangent `along`,
95    /// signed: positive where the surface bends towards its normal
96    /// `normalize(Su × Sv)` (see [`Self::normal`]) that way, the inverse of
97    /// the radius of the circle it bends round there, on that side (see
98    /// [`Self::curvature_radius_along`]).
99    pub fn normal_curvature(&self, u: Scalar, v: Scalar, along: &Vector3) -> GeopResult<Scalar> {
100        let (su, sv) = self.derivatives(u, v)?;
101        let [suu, suv, svv] = self.second_derivatives(u, v)?;
102        let n = su.prod_cross(&sv).normalize()?;
103        let (e, f, g) = (su.prod_dot(&su), su.prod_dot(&sv), sv.prod_dot(&sv));
104        let det = e.mul(g).sub(f.mul(f));
105        let (t_u, t_v) = (along.prod_dot(&su), along.prod_dot(&sv));
106        let ctx = |err: GeopError| {
107            err.with_context(format!(
108                "NurbSurface::normal_curvature({:?}): its first fundamental form E = {e:?}, F = {f:?}, G = {g:?} could be singular",
109                (u, v)
110            ))
111        };
112        let a = g.mul(t_u).sub(f.mul(t_v)).div(det).map_err(ctx)?;
113        let b = e.mul(t_v).sub(f.mul(t_u)).div(det).map_err(ctx)?;
114        Ok(suu
115            .prod_dot(&n)
116            .mul(a)
117            .mul(a)
118            .add(Scalar::TWO.mul(suv.prod_dot(&n)).mul(a).mul(b))
119            .add(svv.prod_dot(&n).mul(b).mul(b)))
120    }
121
122    /// The unit normal `normalize(Su × Sv)` at `(u, v)` and how fast it
123    /// turns as the point moves on the surface at the velocity `along`, a
124    /// tangent vector there: `dn = (dN − n (n·dN)) / |N|` for `N = Su ×
125    /// Sv`, `dN = (Suu a + Suv b) × Sv + Su × (Suv a + Svv b)`, `(a, b)`
126    /// the parameters' rates that move the point at `along` (its
127    /// components in `Su`, `Sv`, by the first fundamental form). Over a
128    /// box of `(u, v)` an enclosure of both over the box.
129    pub fn normal_turning(
130        &self,
131        u: Scalar,
132        v: Scalar,
133        along: &Vector3,
134    ) -> GeopResult<(Vector3, Vector3)> {
135        let s = self.partials(u, v, 2)?;
136        let (su, sv, suu, suv, svv) = (s[1][0], s[0][1], s[2][0], s[1][1], s[0][2]);
137        let (e, f, g) = (su.prod_dot(&su), su.prod_dot(&sv), sv.prod_dot(&sv));
138        let det = e.mul(g).sub(f.mul(f));
139        let (t_u, t_v) = (along.prod_dot(&su), along.prod_dot(&sv));
140        let a = g.mul(t_u).sub(f.mul(t_v)).div(det)?;
141        let b = e.mul(t_v).sub(f.mul(t_u)).div(det)?;
142        let big = su.prod_cross(&sv);
143        let length = big.norm();
144        let n = big.prod_scalar(Scalar::ONE.div(length)?);
145        let turn = suu
146            .prod_scalar(a)
147            .add(&suv.prod_scalar(b))
148            .prod_cross(&sv)
149            .add(&su.prod_cross(&suv.prod_scalar(a).add(&svv.prod_scalar(b))));
150        let dn = turn
151            .sub(&n.prod_scalar(n.prod_dot(&turn)))
152            .prod_scalar(Scalar::ONE.div(length)?);
153        Ok((n, dn))
154    }
155
156    /// The mean curvature at `(u, v)`, signed against the normal
157    /// `normalize(Su × Sv)`: half the trace of the second fundamental form
158    /// relative to the first,
159    /// `H = (L G − 2 M F + N E) / (2 (E G − F²))`, with `E, F, G` the
160    /// first fundamental form and `L, M, N` the second (`Suu·n`, `Suv·n`,
161    /// `Svv·n`). Unlike [`Self::curvature_radius`] it is exact for any
162    /// parametrization, orthogonal or not.
163    ///
164    /// An error where `Su × Sv` vanishes (a pole): the surface may well have
165    /// a curvature there, but not one these partials can express.
166    pub fn mean_curvature(&self, u: Scalar, v: Scalar) -> GeopResult<Scalar> {
167        let (su, sv) = self.derivatives(u, v)?;
168        let [suu, suv, svv] = self.second_derivatives(u, v)?;
169        let cross = su.prod_cross(&sv);
170        let n = cross.normalize()?;
171        let (e, f, g) = (su.prod_dot(&su), su.prod_dot(&sv), sv.prod_dot(&sv));
172        let (l, m, nn) = (suu.prod_dot(&n), suv.prod_dot(&n), svv.prod_dot(&n));
173        l.mul(g)
174            .sub(Scalar::TWO.mul(m).mul(f))
175            .add(nn.mul(e))
176            .div(Scalar::TWO.mul(cross.norm_sq()))
177    }
178}
179
180#[cfg(test)]
181mod tests {
182    use crate::nurb_surface::NurbSurface;
183    use crate::test_fixtures::hom;
184    use geop_core_math::scalars::Scalar;
185
186    /// A flat bilinear patch has zero curvature everywhere: no radius limit.
187    #[test]
188    fn flat_patch_has_no_curvature_radius() {
189        let f = Scalar::from_f64;
190        let s = NurbSurface::try_new(
191            1,
192            1,
193            vec![
194                hom(0., 0., 0., 1.),
195                hom(0., 1., 0., 1.),
196                hom(1., 0., 0., 1.),
197                hom(1., 1., 0., 1.),
198            ],
199            vec![f(0.), f(0.), f(1.), f(1.)],
200            vec![f(0.), f(0.), f(1.), f(1.)],
201        )
202        .unwrap();
203
204        assert!(s.curvature_radius(f(0.5), f(0.5)).unwrap().is_none());
205    }
206
207    /// The saddle `Scalar(u, v) = (u, v, uv)`: its only second partial is the
208    /// mixed one, so it alone gives the mean curvature
209    /// `H = -uv / (1 + u² + v²)^(3/2)`.
210    #[test]
211    fn saddle_mean_curvature() {
212        let f = Scalar::from_f64;
213        let s = NurbSurface::try_new(
214            1,
215            1,
216            vec![
217                hom(0., 0., 0., 1.),
218                hom(0., 1., 0., 1.),
219                hom(1., 0., 0., 1.),
220                hom(1., 1., 1., 1.),
221            ],
222            vec![f(0.), f(0.), f(1.), f(1.)],
223            vec![f(0.), f(0.), f(1.), f(1.)],
224        )
225        .unwrap();
226        let (u, v): (f64, f64) = (0.5, 0.25);
227        let expected = -u * v / (1.0 + u * u + v * v).powf(1.5);
228        let h = s.mean_curvature(f(u), f(v)).unwrap();
229        assert!(h.could_be_equal(f(expected)), "{h:?} vs {expected}");
230    }
231
232    /// A rational quarter cylinder of radius 2 has `|H| = 1 / 4` everywhere,
233    /// where the weights make its parametrization uneven too.
234    #[test]
235    fn cylinder_mean_curvature() {
236        let f = Scalar::from_f64;
237        let (r, w) = (2.0, std::f64::consts::FRAC_1_SQRT_2);
238        let mut cps = Vec::new();
239        for (x, y, wi) in [(1., 0., 1.), (1., 1., w), (0., 1., 1.)] {
240            for z in [0., 1.] {
241                cps.push(hom(x * r * wi, y * r * wi, z * wi, wi));
242            }
243        }
244        let s = NurbSurface::try_new(
245            2,
246            1,
247            cps,
248            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
249            vec![f(0.), f(0.), f(1.), f(1.)],
250        )
251        .unwrap();
252        for (u, v) in [(0.1, 0.2), (0.5, 0.5), (0.8, 0.9)] {
253            let h = s.mean_curvature(f(u), f(v)).unwrap();
254            assert!(h.abs().could_be_equal(f(0.25)), "{h:?} at ({u}, {v})");
255        }
256    }
257
258    /// A rational quarter cylinder of radius 2 bends round its axis with
259    /// radius 2 and runs straight along it — and every direction between
260    /// bends with a radius of `2 / sin² θ` — where the tightest bend, by
261    /// [`NurbSurface::curvature_radius`], is the one round it.
262    #[test]
263    fn a_cylinder_bends_as_its_direction_says() {
264        use geop_core_math::vector::Vector3;
265        let f = Scalar::from_f64;
266        let (r, w) = (2.0, std::f64::consts::FRAC_1_SQRT_2);
267        let mut cps = Vec::new();
268        for (x, y, wi) in [(1., 0., 1.), (1., 1., w), (0., 1., 1.)] {
269            for z in [0., 1.] {
270                cps.push(hom(x * r * wi, y * r * wi, z * wi, wi));
271            }
272        }
273        let s = NurbSurface::try_new(
274            2,
275            1,
276            cps,
277            vec![f(0.), f(0.), f(0.), f(1.), f(1.), f(1.)],
278            vec![f(0.), f(0.), f(1.), f(1.)],
279        )
280        .unwrap();
281        let (u, v) = (f(0.5), f(0.5));
282        // At u = 1/2 the surface point is at 45 degrees: round is along
283        // (-sin, cos, 0), the axis along z.
284        let h = std::f64::consts::FRAC_1_SQRT_2;
285        let direction = |x: f64, y: f64, z: f64| Vector3::from_array([f(x), f(y), f(z)]);
286        let axis = s
287            .curvature_radius_along(u, v, &direction(0.0, 0.0, 1.0))
288            .unwrap();
289        assert!(axis.is_none(), "{axis:?}");
290        let round = s
291            .curvature_radius_along(u, v, &direction(-h, h, 0.0))
292            .unwrap()
293            .expect("a bend");
294        assert!(round.could_be_equal(f(r)), "{round:?}");
295        // Halfway between: sin² of 45 degrees is one half.
296        let between = s
297            .curvature_radius_along(u, v, &direction(-0.5, 0.5, h))
298            .unwrap()
299            .expect("a bend");
300        assert!(between.could_be_equal(f(2.0 * r)), "{between:?}");
301        let tightest = s.curvature_radius(u, v).unwrap().expect("a bend");
302        assert!(tightest.could_be_equal(f(r)), "{tightest:?}");
303    }
304}
geop_core_geometry/nurb_surface/curvature.rs

geop_core_geometry/nurb_surface/
curvature.rs
