geop_core_sketch/space/
curves.rs

1//! A 3-D sketch's curves as NURBS in space, and its curves joined into
2//! chains — the paths and rails built on it.
3
4use std::collections::{BTreeMap, BTreeSet};
5
6use geop_core_geometry::nurb_curve::{Handedness, Helix, NurbCurve, NurbCurve3D};
7use geop_core_math::{
8    geop_error::{GeopError, GeopResult, WithContext},
9    primitives::CoordinateSystem,
10    scalars::Scalar,
11    solvers::walk::{Walked, walk_by},
12    vector::Vector3,
13    with_context,
14};
15
16use super::{
17    Constraint3d, Curve3d, CurveKind3d, Enclosure3d, End, Sketch3d,
18    geometry::{Arc3, Helix3},
19};
20use crate::{CurveId, PointId, ProfileEdge, ProfileJoint};
21
22/// One NURBS piece of a chain, on `[0, 1]`, named as a planar sketch's
23/// pieces are (see [`crate::profile::Piece`]).
24pub type Piece3d = crate::profile::Piece<NurbCurve3D>;
25
26/// Curves of a 3-D sketch joined end to end: an open chain, or a closed one
27/// whose last curve ends where the first starts. Each curve runs forwards
28/// or, `reversed`, backwards along the chain.
29#[derive(Clone, Debug, PartialEq, Eq)]
30pub struct Chain3d {
31    pub edges: Vec<ProfileEdge>,
32    pub closed: bool,
33}
34
35/// The rational quadratic of the arc from `p` to `q` turning by twice the
36/// angle with cosine `cos` and sine `sin` about `normal`: its middle control
37/// point where the tangents at its ends meet, `(L/2) tan` out from the
38/// chord, towards the arc — exact for any radius, a nearly straight arc's
39/// too.
40fn conic(
41    p: &Vector3,
42    q: &Vector3,
43    normal: &Vector3,
44    cos: Scalar,
45    sin: Scalar,
46) -> GeopResult<NurbCurve3D> {
47    let (mid, out) = chord_frame(p, q, normal)?;
48    let m = mid.add(&out.prod_scalar(sin.div(cos)?));
49    NurbCurve3D::arc(*p, m, *q, cos)
50}
51
52/// The middle of the chord from `p` to `q`, and half its length towards the
53/// arc about `normal` through them — whichever its sweep: the chord turned a
54/// quarter back about the normal.
55fn chord_frame(p: &Vector3, q: &Vector3, normal: &Vector3) -> GeopResult<(Vector3, Vector3)> {
56    let half = Scalar::ONE.div(Scalar::TWO)?;
57    let chord = q.sub(p);
58    let out = chord
59        .prod_cross(normal)
60        .normalize()?
61        .prod_scalar(chord.try_norm()?.mul(half));
62    Ok((p.add(q).prod_scalar(half), out))
63}
64
65/// How often an arc whose half sweep has the cosine `cos` is halved so that
66/// each piece's middle weight stays positive: not at all up to a quarter
67/// turn, once up to a half turn, twice beyond — one, two or four pieces.
68pub(super) fn halvings(cos: Scalar) -> GeopResult<u32> {
69    Ok(
70        if !cos.definitely_less(Scalar::ONE.div(Scalar::TWO)?.sqrt()?) {
71            0
72        } else if !cos.definitely_less(Scalar::ZERO) {
73            1
74        } else {
75            2
76        },
77    )
78}
79
80/// The arc from `from` to `to` about `normal`, `cos` and `sin` of half its
81/// sweep, as rational quadratics: halved `halvings` times, each piece
82/// bisected at its middle, `(L/2) tan(quarter)` out from its chord.
83pub(super) fn arc_pieces(
84    from: &Vector3,
85    to: &Vector3,
86    normal: &Vector3,
87    (cos, sin): (Scalar, Scalar),
88    halvings: u32,
89) -> GeopResult<Vec<NurbCurve3D>> {
90    let mut ends = vec![*from, *to];
91    let (mut cos, mut sin) = (cos, sin);
92    for _ in 0..halvings {
93        let tan_quarter = sin.div(Scalar::ONE.add(cos))?;
94        let mut halved = vec![ends[0]];
95        for w in ends.windows(2) {
96            let (mid, out) = chord_frame(&w[0], &w[1], normal)?;
97            halved.push(mid.add(&out.prod_scalar(tan_quarter)));
98            halved.push(w[1]);
99        }
100        ends = halved;
101        let half_cos = Scalar::ONE.add(cos).div(Scalar::TWO)?.sqrt()?;
102        sin = sin.div(Scalar::TWO.mul(half_cos))?;
103        cos = half_cos;
104    }
105    ends.windows(2)
106        .map(|w| conic(&w[0], &w[1], normal, cos, sin))
107        .collect()
108}
109
110impl Sketch3d {
111    /// The curve `curve` as NURBS pieces from its start to its end, each on
112    /// `[0, 1]`, built from `geometry`. How many pieces — an arc has one per
113    /// quarter turn, or two or four, a rational quadratic's middle weight
114    /// staying positive — is decided on the sketch as drawn, which is
115    /// design data, so the pieces' names do not depend on how precisely the
116    /// geometry is known.
117    pub fn curve_nurbs(
118        &self,
119        curve: CurveId,
120        geometry: &Enclosure3d,
121    ) -> GeopResult<Vec<NurbCurve3D>> {
122        let ctx = with_context!("3-D sketch curve {curve} as NURBS");
123        let at = |p: &PointId| geometry.points[p];
124        match &self.curve(curve).with_context(ctx)?.kind {
125            CurveKind3d::Line { start, end } => Ok(vec![
126                NurbCurve3D::line(at(start), at(end)).with_context(ctx)?,
127            ]),
128            CurveKind3d::Arc {
129                start,
130                through,
131                end,
132            } => {
133                let drawn = Arc3 {
134                    s: self.points[start].at,
135                    m: self.points[through].at,
136                    e: self.points[end].at,
137                };
138                let arc = Arc3 {
139                    s: at(start),
140                    m: at(through),
141                    e: at(end),
142                };
143                arc_pieces(
144                    &arc.s,
145                    &arc.e,
146                    &arc.normal().with_context(ctx)?,
147                    arc.half_sweep().with_context(ctx)?,
148                    halvings(drawn.half_sweep().with_context(ctx)?.0)?,
149                )
150                .with_context(ctx)
151            }
152            CurveKind3d::Spline { points } => {
153                let through: Vec<Vector3> = points.iter().map(at).collect();
154                let start = self
155                    .end_tangent(curve, End::Start, geometry)
156                    .with_context(ctx)?;
157                let end = self
158                    .end_tangent(curve, End::End, geometry)
159                    .with_context(ctx)?;
160                Ok(vec![
161                    NurbCurve::cubic_spline(&through, start, end).with_context(ctx)?,
162                ])
163            }
164            CurveKind3d::Helix {
165                start,
166                base,
167                top,
168                turns,
169                left,
170                ..
171            } => {
172                let helix = Helix3 {
173                    start: at(start),
174                    base: at(base),
175                    top: at(top),
176                    turns: *turns,
177                    left: *left,
178                };
179                let (w, out) = helix.axes().with_context(ctx)?;
180                let radius = out.try_norm().with_context(ctx)?;
181                let u = out.normalize().with_context(ctx)?;
182                let foot = helix.start.sub(&out);
183                let frame =
184                    CoordinateSystem::try_new(foot, u, w.prod_cross(&u), w).with_context(ctx)?;
185                let rise = helix.top.sub(&helix.base).try_norm().with_context(ctx)?;
186                let handedness = if *left {
187                    Handedness::Left
188                } else {
189                    Handedness::Right
190                };
191                let motion = Helix::new(
192                    rise.div(*turns).with_context(ctx)?,
193                    turns.to_f64(),
194                    handedness,
195                );
196                Ok(vec![
197                    NurbCurve3D::helix(&frame, radius, &motion).with_context(ctx)?,
198                ])
199            }
200            CurveKind3d::Reference => Err(GeopError::new(format!(
201                "curve {curve} is a reference to an edge of the part, not drawn in the sketch"
202            ))),
203        }
204    }
205
206    /// The direction the spline `curve` leaves (`end` its start) or arrives
207    /// (its end) along, if it is given one (see [`Sketch3d::adopts`]): a
208    /// [`Constraint3d::TangentTo`]'s, the way the spline runs there; or the
209    /// direction of the curve it goes on smoothly from or to.
210    fn end_tangent(
211        &self,
212        curve: CurveId,
213        end: End,
214        geometry: &Enclosure3d,
215    ) -> GeopResult<Option<Vector3>> {
216        let CurveKind3d::Spline { points } = &self.curve(curve)?.kind else {
217            unreachable!("only a spline's ends take a direction");
218        };
219        for c in self.constraints.values() {
220            match *c {
221                Constraint3d::TangentTo {
222                    curve: k,
223                    end: e,
224                    direction,
225                } if k == curve && e == end => {
226                    // Along the direction either way: the way the spline
227                    // runs there, from its end point to the next, as drawn.
228                    let (from, to) = match end {
229                        End::Start => (points[0], points[1]),
230                        End::End => (points[points.len() - 2], points[points.len() - 1]),
231                    };
232                    let run = self.points[&to].at.sub(&self.points[&from].at);
233                    let sign = if run.prod_dot(&direction).to_f64() < 0.0 {
234                        Scalar::ONE.neg()
235                    } else {
236                        Scalar::ONE
237                    };
238                    return Ok(Some(direction.prod_scalar(sign)));
239                }
240                Constraint3d::Tangent { a, b } if self.adopts(a, b)? == Some((curve, end)) => {
241                    let other = if a == curve { b } else { a };
242                    let (mine, theirs) = if a == curve {
243                        self.shared_end(a, b)?.expect("validated")
244                    } else {
245                        let (ea, eb) = self.shared_end(a, b)?.expect("validated");
246                        (eb, ea)
247                    };
248                    let tangent = self.tangent_at(other, theirs, geometry)?;
249                    // Their way through the joint is mine where one of us
250                    // arrives and the other leaves.
251                    return Ok(Some(if mine == theirs {
252                        tangent.neg()
253                    } else {
254                        tangent
255                    }));
256                }
257                _ => {}
258            }
259        }
260        Ok(None)
261    }
262
263    /// The tangent of the drawn curve `curve` at its end `end`, in its own
264    /// direction.
265    fn tangent_at(&self, curve: CurveId, end: End, geometry: &Enclosure3d) -> GeopResult<Vector3> {
266        let pieces = self.curve_nurbs(curve, geometry)?;
267        match end {
268            End::Start => pieces[0].tangent(Scalar::ZERO),
269            End::End => pieces[pieces.len() - 1].tangent(Scalar::ONE),
270        }
271    }
272
273    /// Points along the curve `curve` as drawn, from its start to its end —
274    /// dense enough to draw it.
275    pub fn curve_polyline(&self, curve: CurveId) -> GeopResult<Vec<Vector3>> {
276        let pieces = self.curve_nurbs(curve, &Enclosure3d::as_drawn(self))?;
277        crate::profile::polyline(&pieces, NurbCurve3D::evaluate)
278    }
279
280    /// The curves that become edges, oldest first: those drawn, but for
281    /// construction curves.
282    pub fn drawn_curves(&self) -> impl Iterator<Item = (&CurveId, &Curve3d)> {
283        self.curves
284            .iter()
285            .filter(|(_, c)| c.is_drawn() && !c.construction)
286    }
287
288    /// The drawn curves — not construction, not references — joined into
289    /// chains where they share a point (see [`Sketch3d::point_classes`]).
290    ///
291    /// Each chain runs the way its oldest curve (the lowest id) does; a
292    /// closed one starts at that curve, an open one at its end before it.
293    /// Chains come oldest first. Fails where three curves or more meet at a
294    /// point — a branch is no chain — or a single curve closes on itself.
295    pub fn chains(&self) -> GeopResult<Vec<Chain3d>> {
296        let class = self.point_classes();
297        let drawn: Vec<(CurveId, PointId, PointId)> = self
298            .drawn_curves()
299            .filter_map(|(&id, c)| {
300                let (s, e) = c.endpoints()?;
301                Some((id, class[&s], class[&e]))
302            })
303            .collect();
304        let mut at: BTreeMap<PointId, Vec<CurveId>> = BTreeMap::new();
305        for &(id, s, e) in &drawn {
306            if s == e {
307                return Err(GeopError::new(format!(
308                    "curve {id} ends where it starts: a closed chain needs at least two curves"
309                )));
310            }
311            at.entry(s).or_default().push(id);
312            at.entry(e).or_default().push(id);
313        }
314        if let Some((p, curves)) = at.iter().find(|(_, c)| c.len() > 2) {
315            return Err(GeopError::new(format!(
316                "curves {curves:?} meet at point {p}: a chain does not branch"
317            )));
318        }
319        let ends: BTreeMap<CurveId, (PointId, PointId)> =
320            drawn.iter().map(|&(id, s, e)| (id, (s, e))).collect();
321        // The other curve at the point `p` of the curve `c`.
322        let next = |c: CurveId, p: PointId| at[&p].iter().copied().find(|&k| k != c);
323        // Where a chain is not a path: no curve meets more than one other
324        // at a point, so it cannot.
325        let runs_into_itself =
326            |k: CurveId| GeopError::new(format!("the chain through curve {k} runs into itself"));
327        let mut used = BTreeSet::new();
328        let mut chains = Vec::new();
329        for &(first, _, _) in &drawn {
330            if used.contains(&first) {
331                continue;
332            }
333            // Back from the oldest curve's start to where the chain begins:
334            // each curve, and the point it is left at.
335            let back = walk_by(
336                (first, ends[&first].0),
337                |&(c, _)| c,
338                |&(c, p)| {
339                    Ok::<_, GeopError>(next(c, p).map(|k| {
340                        let (s, e) = ends[&k];
341                        (k, if s == p { e } else { s })
342                    }))
343                },
344            )?;
345            let closed = match back.end {
346                Walked::Closed => true,
347                Walked::Open => false,
348                Walked::Revisited((k, _)) => return Err(runs_into_itself(k)),
349            };
350            // Then on from there: each curve, and the point it is entered at.
351            let start = if closed {
352                (first, ends[&first].0)
353            } else {
354                *back.nodes.last().expect("never empty")
355            };
356            let on = walk_by(
357                start,
358                |&(c, _)| c,
359                |&(c, from)| {
360                    let (s, e) = ends[&c];
361                    let to = if s != from { s } else { e };
362                    Ok::<_, GeopError>(next(c, to).map(|k| (k, to)))
363                },
364            )?;
365            if let Walked::Revisited((k, _)) = on.end {
366                return Err(runs_into_itself(k));
367            }
368            let edges: Vec<ProfileEdge> = (on.nodes.iter())
369                .map(|&(curve, from)| ProfileEdge {
370                    curve,
371                    reversed: ends[&curve].0 != from,
372                    contour: 0,
373                })
374                .collect();
375            used.extend(edges.iter().map(|edge| edge.curve));
376            chains.push(Chain3d { edges, closed });
377        }
378        Ok(chains)
379    }
380}
381
382impl Chain3d {
383    /// The chain as NURBS pieces, each on `[0, 1]` and starting exactly
384    /// where the one before ends — the sketch's points at `geometry` —
385    /// named after the sketch's curves and points (see [`Piece3d`]).
386    pub fn to_nurbs(&self, sketch: &Sketch3d, geometry: &Enclosure3d) -> GeopResult<Vec<Piece3d>> {
387        let mut out = Vec::new();
388        for edge in &self.edges {
389            let curve = edge.curve;
390            let pieces = sketch.curve_nurbs(curve, geometry)?;
391            let (first, last) = sketch
392                .curve(curve)?
393                .endpoints()
394                .ok_or_else(|| GeopError::new(format!("curve {curve} of a chain has no ends")))?;
395            let n = pieces.len();
396            let joint = |j: usize| match j {
397                0 => ProfileJoint::Point(first),
398                j if j == n => ProfileJoint::Point(last),
399                index => ProfileJoint::Split { curve, index },
400            };
401            let pieces = pieces.into_iter().enumerate().map(|(index, c)| Piece3d {
402                curve: c,
403                source: curve,
404                index,
405                start: joint(index),
406                end: joint(index + 1),
407            });
408            if edge.reversed {
409                out.extend(pieces.rev().map(|p| Piece3d {
410                    curve: p.curve.reverse(),
411                    start: p.end,
412                    end: p.start,
413                    ..p
414                }));
415            } else {
416                out.extend(pieces);
417            }
418        }
419        Ok(out)
420    }
421}
geop_core_sketch/space/curves.rs

geop_core_sketch/space/
curves.rs
