axiolid_nurbs/
surface_analysis.rs

1//! Differential properties of regular parametric surfaces.
2
3use axiolid_contracts::{GeomError, GeomResult};
4use axiolid_core::{Point3, Scalar, Tolerance, Vec3};
5use axiolid_evaluate::surface::jet;
6use axiolid_surface::Surface;
7
8/// Coefficients `(e, f, g)` of a quadratic fundamental form.
9#[derive(Debug, Clone, Copy, PartialEq)]
10pub struct FundamentalForm {
11    /// First diagonal coefficient (`E` for the first form, `e` for the second).
12    pub e: Scalar,
13    /// Mixed coefficient (`F` for the first form, `f` for the second).
14    pub f: Scalar,
15    /// Second diagonal coefficient (`G` for the first form, `g` for the second).
16    pub g: Scalar,
17}
18
19/// First/second-order differential properties of an oriented surface.
20#[derive(Debug, Clone, Copy, PartialEq)]
21pub struct SurfaceDifferential {
22    /// Evaluated surface point.
23    pub point: Point3,
24    /// Unit normal from `du × dv`.
25    pub unit_normal: Vec3,
26    /// First fundamental form `(E, F, G)`.
27    pub first: FundamentalForm,
28    /// Second fundamental form `(e, f, g)` for `unit_normal`.
29    pub second: FundamentalForm,
30    /// Gaussian curvature, invariant under normal reversal.
31    pub gaussian_curvature: Scalar,
32    /// Mean curvature for the reported normal orientation.
33    pub mean_curvature: Scalar,
34    /// Principal curvatures in ascending order for the reported orientation.
35    pub principal_curvatures: [Scalar; 2],
36}
37
38/// Analyze a regular oriented surface at `(u, v)`.
39pub fn analyze_surface(
40    surface: &Surface,
41    u: Scalar,
42    v: Scalar,
43    tolerance: Tolerance,
44) -> GeomResult<SurfaceDifferential> {
45    let j = jet(surface, u, v)?;
46    let cross = j.du.cross(j.dv);
47    let area = cross.length();
48    let area_floor = tolerance.linear() * tolerance.linear();
49    if !area.is_finite() || area <= area_floor {
50        return Err(GeomError::Degenerate(format!(
51            "surface differential area {area} does not exceed the tolerance area"
52        )));
53    }
54    let n = cross / area;
55    let first = FundamentalForm {
56        e: j.du.dot(j.du),
57        f: j.du.dot(j.dv),
58        g: j.dv.dot(j.dv),
59    };
60    let second = FundamentalForm {
61        e: n.dot(j.duu),
62        f: n.dot(j.duv),
63        g: n.dot(j.dvv),
64    };
65    let determinant = first.e * first.g - first.f * first.f;
66    if !determinant.is_finite() || determinant <= area_floor * area_floor {
67        return Err(GeomError::Degenerate(
68            "surface first fundamental form is singular".to_owned(),
69        ));
70    }
71    let gaussian = (second.e * second.g - second.f * second.f) / determinant;
72    let mean =
73        (second.e * first.g - 2.0 * second.f * first.f + second.g * first.e) / (2.0 * determinant);
74    let raw_discriminant = mean * mean - gaussian;
75    let scale = mean.abs().max(gaussian.abs().sqrt()).max(1.0);
76    let floor = -64.0 * Scalar::EPSILON * scale * scale;
77    if raw_discriminant < floor || !raw_discriminant.is_finite() {
78        return Err(GeomError::Degenerate(
79            "surface principal-curvature discriminant is invalid".to_owned(),
80        ));
81    }
82    let root = raw_discriminant.max(0.0).sqrt();
83    let result = SurfaceDifferential {
84        point: j.point,
85        unit_normal: n,
86        first,
87        second,
88        gaussian_curvature: gaussian,
89        mean_curvature: mean,
90        principal_curvatures: [mean - root, mean + root],
91    };
92    if [
93        gaussian,
94        mean,
95        result.principal_curvatures[0],
96        result.principal_curvatures[1],
97    ]
98    .iter()
99    .all(|x| x.is_finite())
100    {
101        Ok(result)
102    } else {
103        Err(GeomError::Degenerate(
104            "surface curvature is non-finite".to_owned(),
105        ))
106    }
107}
axiolid_nurbs/surface_analysis.rs

axiolid_nurbs/
surface_analysis.rs
