ogeom_heal/
remodel.rs

1//! Remodelling: a shape's geometry restated wholesale, every trim
2//! re-derived against what replaced its surface.
3//!
4//! [`restrict_degree`] brings every spline curve and surface under a degree
5//! limit at a stated tolerance, for a format that caps it.
6//! [`swept_to_elementary`] names the swept surfaces that are planes, drums,
7//! cones, balls or tori as what they are: an extruded line or circle, a
8//! revolved line or circle.
9
10use ogeom_algo::{Built, normals_oppose, project_on_surface, restate_geometry};
11use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
12use ogeom_geom::{
13    ConeSurface, Curve, CylinderSurface, PlaneSurface, SphereSurface, Surface as _,
14    SurfaceGeometry, TorusSurface,
15};
16use ogeom_math::{Axis, Cone, Cylinder, Direction, Frame, Plane, Point, Sphere, Torus, Vector};
17use ogeom_topo::{Model, Shape};
18
19/// `shape` with every spline curve and surface of degree above
20/// `max_degree` refitted at that degree, each within `tolerance` of what
21/// it replaces, and every trim re-derived. Analytic geometry has no degree
22/// to restrict and stays as it is.
23///
24/// # Errors
25///
26/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
27/// `max_degree` is zero or `tolerance` is not a distance;
28/// [`OgeomError::NotDone`](ogeom_core::OgeomError::NotDone) if a fit cannot
29/// reach `tolerance`; as [`ogeom_algo::restate_geometry`].
30pub fn restrict_degree(
31    model: &mut Model,
32    shape: &Shape,
33    max_degree: usize,
34    tolerance: f64,
35    tol: Tolerances,
36) -> OgeomResult<Built> {
37    if max_degree == 0 {
38        ogeom_bail!(Construction, "a degree limit of zero leaves nothing");
39    }
40    if !(tolerance.is_finite() && tolerance > 0.0) {
41        ogeom_bail!(Construction, "a tolerance of {tolerance} is not a distance");
42    }
43    let surface = |s: &SurfaceGeometry| -> OgeomResult<Option<(SurfaceGeometry, bool)>> {
44        let SurfaceGeometry::BSpline(b) = s else {
45            return Ok(None);
46        };
47        if b.u_knots().degree() <= max_degree && b.v_knots().degree() <= max_degree {
48            return Ok(None);
49        }
50        let fitted = b.restricted_to_degree(max_degree, tolerance, tol)?;
51        if !fitted.met {
52            ogeom_bail!(
53                NotDone,
54                "a patch fitted at degree {max_degree} stays {} away",
55                fitted.error
56            );
57        }
58        let restated = SurfaceGeometry::BSpline(fitted.curve);
59        let flipped = normals_oppose(s, &restated, tol)?;
60        Ok(Some((restated, flipped)))
61    };
62    let curve = |c: &Curve, range: (f64, f64)| -> OgeomResult<Option<(Curve, (f64, f64))>> {
63        let Curve::BSpline(b) = c else {
64            return Ok(None);
65        };
66        if b.degree() <= max_degree {
67            return Ok(None);
68        }
69        // The edge's piece, open: a ring's spline has no ends to fit
70        // between, and its piece does.
71        let piece = b.segment(range, tol)?;
72        let fitted = piece.restricted_to_degree(max_degree, tolerance, tol)?;
73        if !fitted.met {
74            ogeom_bail!(
75                NotDone,
76                "a curve fitted at degree {max_degree} stays {} away",
77                fitted.error
78            );
79        }
80        let restated = Curve::BSpline(fitted.curve);
81        let domain = ogeom_geom::Curve3d::domain(&restated);
82        Ok(Some((restated, domain)))
83    };
84    restate_geometry(model, shape, &surface, &curve, tol)
85}
86
87/// `shape` with every surface of extrusion or revolution that is a plane,
88/// drum, cone, ball or torus restated as it: a line or circle extruded, a
89/// line or circle revolved in a plane through the axis. Every trim is
90/// re-derived; the surfaces are the same point sets, so nothing moves.
91///
92/// # Errors
93///
94/// As [`ogeom_algo::restate_geometry`].
95pub fn swept_to_elementary(
96    model: &mut Model,
97    shape: &Shape,
98    tol: Tolerances,
99) -> OgeomResult<Built> {
100    let surface = |s: &SurfaceGeometry| -> OgeomResult<Option<(SurfaceGeometry, bool)>> {
101        let Some(candidate) = elementary(s, tol)? else {
102            return Ok(None);
103        };
104        // Held to every sample: a wrong yes is a solid with the wrong
105        // surface under every later operation.
106        for p in interior_samples(s, 5, tol)? {
107            if project_on_surface(&candidate, p, 16, tol)?.distance > tol.confusion() {
108                return Ok(None);
109            }
110        }
111        let flipped = normals_oppose(s, &candidate, tol)?;
112        Ok(Some((candidate, flipped)))
113    };
114    let curve = |_: &Curve, _: (f64, f64)| -> OgeomResult<Option<(Curve, (f64, f64))>> { Ok(None) };
115    restate_geometry(model, shape, &surface, &curve, tol)
116}
117
118/// The elementary surface a swept one proposes, unverified.
119fn elementary(s: &SurfaceGeometry, tol: Tolerances) -> OgeomResult<Option<SurfaceGeometry>> {
120    let samples = interior_samples(s, 5, tol)?;
121    let unwrap = |c: &Curve| -> Curve {
122        match c {
123            Curve::Trimmed(t) => t.basis().clone(),
124            other => other.clone(),
125        }
126    };
127    let found = match s {
128        SurfaceGeometry::Extrusion(e) => match unwrap(e.curve()) {
129            Curve::Line(l) => {
130                let axis = l.axis();
131                let normal = axis.direction.vector().cross(e.direction().vector());
132                let Ok(z) = Direction::new(normal, tol) else {
133                    return Ok(None);
134                };
135                let frame = Frame::new(axis.location, z, axis.direction, tol)?;
136                Some(plane_over(frame, &samples)?)
137            }
138            Curve::Circle(c) => {
139                let circle = c.circle();
140                let frame = circle.frame();
141                if frame.z().vector().cross(e.direction().vector()).magnitude() > tol.angular() {
142                    return Ok(None);
143                }
144                let axis = Frame::new(frame.origin(), e.direction(), frame.x(), tol)?;
145                let heights = heights(&axis, &samples);
146                Some(
147                    CylinderSurface::new(Cylinder::new(axis, circle.radius(), tol)?, heights)?
148                        .into(),
149                )
150            }
151            _ => None,
152        },
153        SurfaceGeometry::Revolution(r) => {
154            let axis = r.axis();
155            match unwrap(r.curve()) {
156                Curve::Line(l) => revolved_line(axis, l.axis(), &samples, tol)?,
157                Curve::Circle(c) => revolved_circle(axis, c.circle(), tol)?,
158                _ => None,
159            }
160        }
161        _ => None,
162    };
163    Ok(found)
164}
165
166/// A line revolved in a plane through the axis: a drum where it runs with
167/// the axis, a disc where it runs square to it, a cone otherwise.
168fn revolved_line(
169    axis: Axis,
170    line: Axis,
171    samples: &[Point],
172    tol: Tolerances,
173) -> OgeomResult<Option<SurfaceGeometry>> {
174    let z = axis.direction.vector();
175    let a = line.location;
176    let b = a + line.direction.vector();
177    let height = |p: Point| (p - axis.location).dot(z);
178    let radius = |p: Point| p.distance(axis.project(p));
179    // The profile's own start side names the chart's `u = 0`.
180    let Some(off) = [a, b]
181        .into_iter()
182        .map(|p| p - axis.project(p))
183        .find(|v| v.magnitude() > tol.confusion())
184    else {
185        return Ok(None);
186    };
187    let x = Direction::new(off, tol)?;
188    let (dh, dr) = (height(b) - height(a), radius(b) - radius(a));
189    if dh.abs() <= tol.confusion() {
190        let frame = Frame::new(axis.project(a), axis.direction, x, tol)?;
191        return Ok(Some(plane_over(frame, samples)?));
192    }
193    let frame = Frame::new(axis.project(a), axis.direction, x, tol)?;
194    let heights = heights(&frame, samples);
195    if dr.abs() <= tol.confusion() {
196        return Ok(Some(
197            CylinderSurface::new(Cylinder::new(frame, radius(a), tol)?, heights)?.into(),
198        ));
199    }
200    let half_angle = (dr / dh).atan();
201    let Ok(cone) = Cone::new(frame, radius(a), half_angle, tol) else {
202        return Ok(None);
203    };
204    Ok(Some(ConeSurface::new(cone, heights)?.into()))
205}
206
207/// A circle revolved in a plane through the axis: a ball where its centre
208/// is on the axis, a torus otherwise.
209fn revolved_circle(
210    axis: Axis,
211    circle: ogeom_math::Circle,
212    tol: Tolerances,
213) -> OgeomResult<Option<SurfaceGeometry>> {
214    let centre = circle.centre();
215    let foot = axis.project(centre);
216    let off = centre - foot;
217    if off.magnitude() <= tol.confusion() {
218        let x = Direction::new(circle.frame().x().vector(), tol)?;
219        let x = if x.vector().cross(axis.direction.vector()).magnitude() > tol.angular() {
220            x
221        } else {
222            Direction::new(circle.frame().y().vector(), tol)?
223        };
224        let side = x.vector() - axis.direction.vector() * x.vector().dot(axis.direction.vector());
225        let frame = Frame::new(foot, axis.direction, Direction::new(side, tol)?, tol)?;
226        return Ok(Some(
227            SphereSurface::new(Sphere::new(frame, circle.radius(), tol)?).into(),
228        ));
229    }
230    let frame = Frame::new(foot, axis.direction, Direction::new(off, tol)?, tol)?;
231    let Ok(torus) = Torus::new(frame, off.magnitude(), circle.radius(), tol) else {
232        return Ok(None);
233    };
234    Ok(Some(TorusSurface::new(torus).into()))
235}
236
237/// A plane in `frame`, bounded round the samples.
238fn plane_over(frame: Frame, samples: &[Point]) -> OgeomResult<SurfaceGeometry> {
239    let along = |p: Point, d: Vector| (p - frame.origin()).dot(d);
240    let range = |d: Vector| {
241        let (lo, hi) = samples
242            .iter()
243            .fold((f64::INFINITY, f64::NEG_INFINITY), |(l, h), p| {
244                (l.min(along(*p, d)), h.max(along(*p, d)))
245            });
246        let margin = (hi - lo) * 0.5 + 1.0;
247        (lo - margin, hi + margin)
248    };
249    Ok(PlaneSurface::over(
250        Plane::new(frame),
251        range(frame.x().vector()),
252        range(frame.y().vector()),
253    )?
254    .into())
255}
256
257/// The heights along a frame's `z` that the samples span, with room.
258fn heights(frame: &Frame, samples: &[Point]) -> (f64, f64) {
259    let z = frame.z().vector();
260    let (lo, hi) = samples
261        .iter()
262        .map(|p| (*p - frame.origin()).dot(z))
263        .fold((f64::INFINITY, f64::NEG_INFINITY), |(l, h), v| {
264            (l.min(v), h.max(v))
265        });
266    let margin = (hi - lo) * 0.5 + 1.0;
267    (lo - margin, hi + margin)
268}
269
270/// Points over the inside of a surface's chart, corners included.
271fn interior_samples(s: &SurfaceGeometry, n: usize, tol: Tolerances) -> OgeomResult<Vec<Point>> {
272    let ((u0, u1), (v0, v1)) = s.domain();
273    let mut out = Vec::with_capacity((n + 1) * (n + 1));
274    for i in 0..=n {
275        for j in 0..=n {
276            #[allow(clippy::cast_precision_loss, reason = "a sample index")]
277            let (fu, fv) = (i as f64 / n as f64, j as f64 / n as f64);
278            out.push(s.point_at(u0 + (u1 - u0) * fu, v0 + (v1 - v0) * fv, tol)?);
279        }
280    }
281    Ok(out)
282}
ogeom_heal/remodel.rs

ogeom_heal/
remodel.rs
