ogeom_geom/
curve2d.rs

1//! Curves in the plane: pcurves.
2//!
3//! These carry a curve through a surface's `(u, v)` parameter space. An edge
4//! holds one per adjacent face (`docs/DATA_MODEL.md` §6), and boolean face
5//! splitting happens entirely in this space: without a pcurve on each face
6//! there is nothing to split *with*.
7//!
8//! A separate type hierarchy from [`crate::curve`] rather than a generic
9//! parameter, because a pcurve is used differently from a spatial curve.
10//! Distance in parameter space is not distance in space; the same parametric
11//! step covers a metre near a cylinder's equator and nothing at all near a
12//! sphere's pole, so a function that treats the two alike is wrong, and
13//! separate types keep that from compiling.
14
15use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
16use ogeom_math::{
17    Axis2, Circle2, Direction2, Ellipse2, KnotVector, Point2, Transform2, Vector2, Weighted,
18    bspline, elementary,
19};
20
21use crate::traits::{Curve2d, CurveKind, Reversible};
22
23const TAU: f64 = core::f64::consts::TAU;
24
25/// A curve in the plane.
26#[derive(Debug, Clone, PartialEq)]
27pub enum PlanarCurve {
28    /// A straight line.
29    Line(Line2d),
30    /// A circle or arc.
31    Circle(Circle2d),
32    /// An ellipse or arc.
33    Ellipse(Ellipse2d),
34    /// A polynomial or rational B-spline.
35    BSpline(BSpline2d),
36    /// Another planar curve restricted to a sub-interval.
37    Trimmed(Box<Trimmed2d>),
38    /// A curve at a constant signed distance along another's left normal.
39    Offset(Box<Offset2d>),
40    /// An affine-plus-trigonometric curve: `c + d·t + a·cos t + b·sin t`:
41    /// the exact chart trace of an oblique analytic section on a periodic
42    /// surface.
43    Trig(Trig2d),
44}
45
46/// A straight line in the plane, parameterized by length.
47#[derive(Debug, Clone, Copy, PartialEq)]
48pub struct Line2d {
49    axis: Axis2,
50    domain: (f64, f64),
51}
52
53/// A circle in the plane, parameterized by angle.
54#[derive(Debug, Clone, Copy, PartialEq)]
55pub struct Circle2d {
56    circle: Circle2,
57    reversed: bool,
58}
59
60/// An ellipse in the plane, parameterized by eccentric angle.
61#[derive(Debug, Clone, Copy, PartialEq)]
62pub struct Ellipse2d {
63    ellipse: Ellipse2,
64    reversed: bool,
65}
66
67/// A B-spline in the plane, polynomial or rational.
68#[derive(Debug, Clone, PartialEq)]
69pub struct BSpline2d {
70    knots: KnotVector,
71    control: Vec<Weighted<Point2>>,
72    rational: bool,
73}
74
75/// Another planar curve restricted to a sub-interval.
76#[derive(Debug, Clone, PartialEq)]
77pub struct Trimmed2d {
78    basis: PlanarCurve,
79    domain: (f64, f64),
80    reversed: bool,
81}
82
83/// Reverse a parameter within `[a, b]`, preserving the interval.
84fn mirror(u: f64, a: f64, b: f64) -> f64 {
85    a + b - u
86}
87
88/// The affine-plus-trigonometric curve `p(t) = c + d·t + a·cos t + b·sin t`.
89///
90/// This is the family chart traces of oblique analytic sections live in: an
91/// oblique plane's ellipse on a cylinder runs linearly in the chart angle
92/// and sinusoidally in height, which no line, conic or spline states
93/// exactly. Closed under similarity transforms and reversal, exact
94/// derivatives to any order.
95#[derive(Debug, Clone, Copy, PartialEq)]
96pub struct Trig2d {
97    c: Point2,
98    d: Vector2,
99    a: Vector2,
100    b: Vector2,
101    domain: (f64, f64),
102    reversed: bool,
103}
104
105impl Trig2d {
106    /// A trig curve over an increasing domain.
107    ///
108    /// # Errors
109    ///
110    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
111    /// any coefficient is not finite or the domain is not increasing.
112    pub fn new(
113        c: Point2,
114        d: Vector2,
115        a: Vector2,
116        b: Vector2,
117        domain: (f64, f64),
118    ) -> OgeomResult<Self> {
119        let finite = c.to_vector().is_finite()
120            && d.is_finite()
121            && a.is_finite()
122            && b.is_finite()
123            && domain.0.is_finite()
124            && domain.1.is_finite();
125        if !finite || domain.0 >= domain.1 {
126            ogeom_bail!(
127                Construction,
128                "a trig curve needs finite coefficients and an increasing domain"
129            );
130        }
131        Ok(Self {
132            c,
133            d,
134            a,
135            b,
136            domain,
137            reversed: false,
138        })
139    }
140
141    /// The constant term.
142    #[must_use]
143    pub const fn constant(&self) -> Point2 {
144        self.c
145    }
146
147    /// The linear coefficient.
148    #[must_use]
149    pub const fn linear(&self) -> Vector2 {
150        self.d
151    }
152
153    /// The cosine coefficient.
154    #[must_use]
155    pub const fn cosine(&self) -> Vector2 {
156        self.a
157    }
158
159    /// The sine coefficient.
160    #[must_use]
161    pub const fn sine(&self) -> Vector2 {
162        self.b
163    }
164
165    /// Whether evaluation runs the domain backwards.
166    #[must_use]
167    pub const fn is_reversed(&self) -> bool {
168        self.reversed
169    }
170
171    fn raw(&self, t: f64) -> Point2 {
172        let (sin, cos) = t.sin_cos();
173        self.c + self.d * t + self.a * cos + self.b * sin
174    }
175}
176
177impl Curve2d for Trig2d {
178    fn domain(&self) -> (f64, f64) {
179        self.domain
180    }
181
182    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point2> {
183        let u = clamp_to_domain(u, self.domain, false, tol)?;
184        let t = if self.reversed {
185            self.domain.0 + self.domain.1 - u
186        } else {
187            u
188        };
189        Ok(self.raw(t))
190    }
191
192    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector2> {
193        Ok(self.derivatives_at(u, 1, tol)?[1])
194    }
195
196    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector2>> {
197        let u = clamp_to_domain(u, self.domain, false, tol)?;
198        let t = if self.reversed {
199            self.domain.0 + self.domain.1 - u
200        } else {
201            u
202        };
203        let (sin, cos) = t.sin_cos();
204        let sign = if self.reversed { -1.0 } else { 1.0 };
205        let mut out = Vec::with_capacity(n + 1);
206        out.push(self.raw(t).to_vector());
207        for order in 1..=n {
208            // The trig part cycles with period four; the linear part
209            // survives only to first order. Odd orders pick up the
210            // reversal sign.
211            let trig = match order % 4 {
212                1 => self.a * -sin + self.b * cos,
213                2 => self.a * -cos + self.b * -sin,
214                3 => self.a * sin + self.b * -cos,
215                _ => self.a * cos + self.b * sin,
216            };
217            let linear = if order == 1 { self.d } else { Vector2::ZERO };
218            let odd = if order % 2 == 1 { sign } else { 1.0 };
219            out.push((trig + linear) * odd);
220        }
221        Ok(out)
222    }
223
224    fn kind(&self) -> CurveKind {
225        CurveKind::Trig
226    }
227
228    fn is_closed(&self, tol: Tolerances) -> bool {
229        self.raw(self.domain.0).distance(self.raw(self.domain.1)) <= tol.confusion()
230    }
231
232    fn is_periodic(&self) -> bool {
233        false
234    }
235}
236
237/// A planar curve displaced a constant signed distance along its basis's
238/// left normal (the tangent turned a quarter left), sharing the basis's
239/// parameterization.
240///
241/// Point and first derivative are exact from the basis's first and second
242/// derivatives; the second would need the basis's third, which the
243/// vocabulary does not carry, so `d2_at` refuses by name. Where the basis's
244/// tangent vanishes the offset direction is undefined and evaluation
245/// refuses.
246#[derive(Debug, Clone, PartialEq)]
247pub struct Offset2d {
248    basis: PlanarCurve,
249    distance: f64,
250}
251
252impl Offset2d {
253    /// Offset `basis` by a signed `distance` along its left normal.
254    ///
255    /// # Errors
256    ///
257    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
258    /// distance is not finite and non-zero.
259    pub fn new(basis: PlanarCurve, distance: f64) -> OgeomResult<Self> {
260        if !distance.is_finite() || distance == 0.0 {
261            ogeom_bail!(
262                Construction,
263                "an offset of {distance} is not a displacement"
264            );
265        }
266        Ok(Self { basis, distance })
267    }
268
269    /// The curve being offset.
270    #[must_use]
271    pub const fn basis(&self) -> &PlanarCurve {
272        &self.basis
273    }
274
275    /// The signed displacement along the left normal.
276    #[must_use]
277    pub const fn distance(&self) -> f64 {
278        self.distance
279    }
280
281    /// The unit left normal and its derivative at `t`, from the basis's
282    /// first two derivatives.
283    fn normal_and_slope(&self, t: f64, tol: Tolerances) -> OgeomResult<(Vector2, Vector2)> {
284        let d = self.basis.derivatives_at(t, 2, tol)?;
285        let rot = |v: Vector2| Vector2::new(-v.y, v.x);
286        let w = rot(d[1]);
287        let m = w.magnitude();
288        if m <= tol.confusion() {
289            ogeom_bail!(
290                Construction,
291                "the basis has no tangent at {t}; the offset direction is undefined"
292            );
293        }
294        let n = w / m;
295        let wp = rot(d[2]);
296        let np = (wp - n * n.dot(wp)) / m;
297        Ok((n, np))
298    }
299}
300
301impl Curve2d for Offset2d {
302    fn domain(&self) -> (f64, f64) {
303        self.basis.domain()
304    }
305
306    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point2> {
307        let base = self.basis.point_at(u, tol)?;
308        let (n, _) = self.normal_and_slope(u, tol)?;
309        Ok(base + n * self.distance)
310    }
311
312    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector2> {
313        let d1 = self.basis.d1_at(u, tol)?;
314        let (_, np) = self.normal_and_slope(u, tol)?;
315        Ok(d1 + np * self.distance)
316    }
317
318    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector2>> {
319        if n >= 2 {
320            ogeom_bail!(
321                Construction,
322                "an offset curve's second derivative needs its basis's third, which the \
323                 vocabulary does not carry"
324            );
325        }
326        let mut out = vec![self.point_at(u, tol)?.to_vector()];
327        if n >= 1 {
328            out.push(self.d1_at(u, tol)?);
329        }
330        Ok(out)
331    }
332
333    fn kind(&self) -> CurveKind {
334        CurveKind::Offset
335    }
336
337    fn is_closed(&self, tol: Tolerances) -> bool {
338        self.basis.is_closed(tol)
339    }
340
341    fn is_periodic(&self) -> bool {
342        self.basis.is_periodic()
343    }
344}
345
346impl Line2d {
347    /// A segment between two distinct points, parameterized by arc length.
348    ///
349    /// # Errors
350    ///
351    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the points
352    /// coincide.
353    pub fn segment(from: Point2, to: Point2, tol: Tolerances) -> OgeomResult<Self> {
354        Ok(Self {
355            axis: Axis2::through(from, to, tol)?,
356            domain: (0.0, from.distance(to)),
357        })
358    }
359
360    /// A line over an explicit parameter range.
361    ///
362    /// # Errors
363    ///
364    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the range
365    /// is empty or non-finite.
366    pub fn over(axis: Axis2, start: f64, end: f64) -> OgeomResult<Self> {
367        if !start.is_finite() || !end.is_finite() || end <= start {
368            ogeom_bail!(Construction, "line range [{start}, {end}] is empty");
369        }
370        Ok(Self {
371            axis,
372            domain: (start, end),
373        })
374    }
375
376    /// The underlying axis.
377    #[must_use]
378    pub const fn axis(&self) -> Axis2 {
379        self.axis
380    }
381}
382
383impl Circle2d {
384    /// A full circle.
385    #[must_use]
386    pub const fn new(circle: Circle2) -> Self {
387        Self {
388            circle,
389            reversed: false,
390        }
391    }
392
393    /// The underlying circle.
394    #[must_use]
395    pub const fn circle(&self) -> Circle2 {
396        self.circle
397    }
398
399    /// Whether the curve runs backwards along its underlying circle.
400    ///
401    /// Part of the curve's state and not derivable from its circle, so
402    /// anything that has to reproduce this curve exactly (the native format
403    /// above all) needs to be able to read it.
404    #[must_use]
405    pub const fn is_reversed(&self) -> bool {
406        self.reversed
407    }
408}
409
410impl Ellipse2d {
411    /// A full ellipse.
412    #[must_use]
413    pub const fn new(ellipse: Ellipse2) -> Self {
414        Self {
415            ellipse,
416            reversed: false,
417        }
418    }
419
420    /// The underlying ellipse.
421    #[must_use]
422    pub const fn ellipse(&self) -> Ellipse2 {
423        self.ellipse
424    }
425
426    /// Whether the curve runs backwards along its underlying ellipse.
427    ///
428    /// Part of the curve's state and not derivable from its ellipse, so
429    /// anything that has to reproduce this curve exactly (the native format
430    /// above all) needs to be able to read it.
431    #[must_use]
432    pub const fn is_reversed(&self) -> bool {
433        self.reversed
434    }
435}
436
437impl BSpline2d {
438    /// A polynomial B-spline.
439    ///
440    /// # Errors
441    ///
442    /// [`OgeomError::Dimension`](ogeom_core::OgeomError::Dimension) on a shape mismatch.
443    pub fn new(knots: KnotVector, control: Vec<Point2>, tol: Tolerances) -> OgeomResult<Self> {
444        let weighted = control
445            .into_iter()
446            .map(|p| Weighted::new(p, 1.0, tol))
447            .collect::<OgeomResult<Vec<_>>>()?;
448        Self::rational(knots, weighted)
449    }
450
451    /// A rational B-spline.
452    ///
453    /// # Errors
454    ///
455    /// [`OgeomError::Dimension`](ogeom_core::OgeomError::Dimension) on a shape mismatch.
456    pub fn rational(knots: KnotVector, control: Vec<Weighted<Point2>>) -> OgeomResult<Self> {
457        if control.len() != knots.control_point_count() {
458            ogeom_bail!(
459                Dimension,
460                "knot vector describes {} control points, got {}",
461                knots.control_point_count(),
462                control.len()
463            );
464        }
465        let first = control[0].weight;
466        let rational = control
467            .iter()
468            .any(|w| (w.weight - first).abs() > 1e-12 * first.abs());
469        Ok(Self {
470            knots,
471            control,
472            rational,
473        })
474    }
475
476    /// The knot vector.
477    #[must_use]
478    pub const fn knots(&self) -> &KnotVector {
479        &self.knots
480    }
481
482    /// The weighted control points.
483    #[must_use]
484    pub fn control_points(&self) -> &[Weighted<Point2>] {
485        &self.control
486    }
487
488    /// Whether the weights differ.
489    #[must_use]
490    pub const fn is_rational(&self) -> bool {
491        self.rational
492    }
493}
494
495impl Trimmed2d {
496    /// Restrict `basis` to `[start, end]`.
497    ///
498    /// # Errors
499    ///
500    /// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if the range is empty or
501    /// leaves the basis curve's domain.
502    pub fn new(basis: PlanarCurve, start: f64, end: f64, tol: Tolerances) -> OgeomResult<Self> {
503        let (a, b) = basis.domain();
504        if !start.is_finite() || !end.is_finite() || end <= start + tol.parametric() {
505            ogeom_bail!(Domain, "trim range [{start}, {end}] is empty");
506        }
507        if !basis.is_periodic() && (start < a - tol.parametric() || end > b + tol.parametric()) {
508            ogeom_bail!(Domain, "trim range [{start}, {end}] leaves [{a}, {b}]");
509        }
510        Ok(Self {
511            basis,
512            domain: (start, end),
513            reversed: false,
514        })
515    }
516
517    /// The curve being trimmed.
518    #[must_use]
519    pub const fn basis(&self) -> &PlanarCurve {
520        &self.basis
521    }
522
523    /// Whether the curve runs backwards along its underlying curve.
524    ///
525    /// Part of the curve's state and not derivable from its basis curve, so
526    /// anything that has to reproduce this curve exactly (the native format
527    /// above all) needs to be able to read it.
528    #[must_use]
529    pub const fn is_reversed(&self) -> bool {
530        self.reversed
531    }
532
533    /// This curve's parameter mapped onto the basis curve's.
534    fn basis_parameter(&self, u: f64, tol: Tolerances) -> OgeomResult<f64> {
535        let u = clamp_to_domain(u, self.domain, false, tol)?;
536        Ok(if self.reversed {
537            mirror(u, self.domain.0, self.domain.1)
538        } else {
539            u
540        })
541    }
542}
543
544/// Bring `u` into `domain`, wrapping if periodic.
545fn clamp_to_domain(
546    u: f64,
547    domain: (f64, f64),
548    periodic: bool,
549    tol: Tolerances,
550) -> OgeomResult<f64> {
551    let (a, b) = domain;
552    if periodic {
553        return Ok(a + (u - a).rem_euclid(b - a));
554    }
555    if !u.is_finite() || u < a - tol.parametric() || u > b + tol.parametric() {
556        ogeom_bail!(Domain, "parameter {u} outside curve domain [{a}, {b}]");
557    }
558    Ok(u.clamp(a, b))
559}
560
561/// Fill a derivative list to `n + 1` entries, padding with zeros.
562fn pad(mut out: Vec<Vector2>, n: usize) -> Vec<Vector2> {
563    out.resize(n.max(out.len().saturating_sub(1)) + 1, Vector2::ZERO);
564    out.truncate(n + 1);
565    out
566}
567
568impl Curve2d for Line2d {
569    fn domain(&self) -> (f64, f64) {
570        self.domain
571    }
572
573    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point2> {
574        Ok(self
575            .axis
576            .point_at(clamp_to_domain(u, self.domain, false, tol)?))
577    }
578
579    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector2> {
580        clamp_to_domain(u, self.domain, false, tol)?;
581        Ok(self.axis.direction.vector())
582    }
583
584    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector2>> {
585        let p = self.point_at(u, tol)?;
586        Ok(pad(vec![p.to_vector(), self.axis.direction.vector()], n))
587    }
588
589    fn kind(&self) -> CurveKind {
590        CurveKind::Line
591    }
592
593    fn is_closed(&self, _tol: Tolerances) -> bool {
594        false
595    }
596
597    fn is_periodic(&self) -> bool {
598        false
599    }
600}
601
602impl Curve2d for Circle2d {
603    fn domain(&self) -> (f64, f64) {
604        (0.0, TAU)
605    }
606
607    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point2> {
608        Ok(Point2::from_vector(self.derivatives_at(u, 0, tol)?[0]))
609    }
610
611    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector2> {
612        Ok(self.derivatives_at(u, 1, tol)?[1])
613    }
614
615    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector2>> {
616        let u = clamp_to_domain(u, self.domain(), true, tol)?;
617        let angle = if self.reversed { -u } else { u };
618        let f = self.circle.frame();
619        let r = self.circle.radius();
620        let (sin, cos) = angle.sin_cos();
621        let (x, y) = (f.x().vector(), f.y().vector());
622        let point = self.circle.centre() + x * (r * cos) + y * (r * sin);
623        // Each order of the reversal picks up a factor of -1, so odd orders flip.
624        let sign = if self.reversed { -1.0 } else { 1.0 };
625        Ok(pad(
626            vec![
627                point.to_vector(),
628                (x * (-r * sin) + y * (r * cos)) * sign,
629                x * (-r * cos) + y * (-r * sin),
630            ],
631            n,
632        ))
633    }
634
635    fn kind(&self) -> CurveKind {
636        CurveKind::Circle
637    }
638
639    fn is_closed(&self, _tol: Tolerances) -> bool {
640        true
641    }
642
643    fn is_periodic(&self) -> bool {
644        true
645    }
646}
647
648impl Curve2d for Ellipse2d {
649    fn domain(&self) -> (f64, f64) {
650        (0.0, TAU)
651    }
652
653    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point2> {
654        Ok(Point2::from_vector(self.derivatives_at(u, 0, tol)?[0]))
655    }
656
657    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector2> {
658        Ok(self.derivatives_at(u, 1, tol)?[1])
659    }
660
661    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector2>> {
662        let u = clamp_to_domain(u, self.domain(), true, tol)?;
663        let angle = if self.reversed { -u } else { u };
664        let f = self.ellipse.frame();
665        let (a, b) = (self.ellipse.major_radius(), self.ellipse.minor_radius());
666        let (sin, cos) = angle.sin_cos();
667        let (x, y) = (f.x().vector(), f.y().vector());
668        let point = self.ellipse.centre() + x * (a * cos) + y * (b * sin);
669        let sign = if self.reversed { -1.0 } else { 1.0 };
670        Ok(pad(
671            vec![
672                point.to_vector(),
673                (x * (-a * sin) + y * (b * cos)) * sign,
674                x * (-a * cos) + y * (-b * sin),
675            ],
676            n,
677        ))
678    }
679
680    fn kind(&self) -> CurveKind {
681        CurveKind::Ellipse
682    }
683
684    fn is_closed(&self, _tol: Tolerances) -> bool {
685        true
686    }
687
688    fn is_periodic(&self) -> bool {
689        true
690    }
691}
692
693impl Curve2d for BSpline2d {
694    fn domain(&self) -> (f64, f64) {
695        self.knots.domain()
696    }
697
698    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point2> {
699        let u = clamp_to_domain(u, self.domain(), false, tol)?;
700        if self.rational {
701            bspline::evaluate_rational(&self.knots, &self.control, u, tol)
702        } else {
703            Ok(bspline::evaluate(&self.knots, &self.control, u, tol)?.point())
704        }
705    }
706
707    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector2> {
708        Ok(self.derivatives_at(u, 1, tol)?[1])
709    }
710
711    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector2>> {
712        let u = clamp_to_domain(u, self.domain(), false, tol)?;
713        let points = bspline::rational_derivatives(&self.knots, &self.control, u, n, tol)?;
714        Ok(points.into_iter().map(Point2::to_vector).collect())
715    }
716
717    fn kind(&self) -> CurveKind {
718        CurveKind::BSpline
719    }
720
721    fn is_closed(&self, tol: Tolerances) -> bool {
722        self.control[0]
723            .point()
724            .is_equal(self.control[self.control.len() - 1].point(), tol)
725    }
726
727    fn is_periodic(&self) -> bool {
728        false
729    }
730}
731
732impl Curve2d for Trimmed2d {
733    fn domain(&self) -> (f64, f64) {
734        self.domain
735    }
736
737    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point2> {
738        self.basis.point_at(self.basis_parameter(u, tol)?, tol)
739    }
740
741    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector2> {
742        let d = self.basis.d1_at(self.basis_parameter(u, tol)?, tol)?;
743        Ok(if self.reversed { -d } else { d })
744    }
745
746    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector2>> {
747        let t = self.basis_parameter(u, tol)?;
748        let mut out = self.basis.derivatives_at(t, n, tol)?;
749        if self.reversed {
750            for (order, d) in out.iter_mut().enumerate() {
751                if order % 2 == 1 {
752                    *d = -*d;
753                }
754            }
755        }
756        Ok(out)
757    }
758
759    fn kind(&self) -> CurveKind {
760        CurveKind::Trimmed
761    }
762
763    fn is_closed(&self, tol: Tolerances) -> bool {
764        match (self.start(tol), self.end(tol)) {
765            (Ok(a), Ok(b)) => a.is_equal(b, tol),
766            _ => false,
767        }
768    }
769
770    fn is_periodic(&self) -> bool {
771        false
772    }
773}
774
775macro_rules! dispatch {
776    ($self:ident, $c:ident => $body:expr) => {
777        match $self {
778            Self::Line($c) => $body,
779            Self::Circle($c) => $body,
780            Self::Ellipse($c) => $body,
781            Self::BSpline($c) => $body,
782            Self::Trimmed($c) => $body,
783            Self::Offset($c) => $body,
784            Self::Trig($c) => $body,
785        }
786    };
787}
788
789impl Curve2d for PlanarCurve {
790    fn domain(&self) -> (f64, f64) {
791        dispatch!(self, c => c.domain())
792    }
793
794    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point2> {
795        dispatch!(self, c => c.point_at(u, tol))
796    }
797
798    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector2> {
799        dispatch!(self, c => c.d1_at(u, tol))
800    }
801
802    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector2>> {
803        dispatch!(self, c => c.derivatives_at(u, n, tol))
804    }
805
806    fn kind(&self) -> CurveKind {
807        dispatch!(self, c => c.kind())
808    }
809
810    fn is_closed(&self, tol: Tolerances) -> bool {
811        dispatch!(self, c => c.is_closed(tol))
812    }
813
814    fn is_periodic(&self) -> bool {
815        dispatch!(self, c => c.is_periodic())
816    }
817}
818
819impl PlanarCurve {
820    /// This curve moved by a planar similarity.
821    ///
822    /// # Errors
823    ///
824    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the result
825    /// would be degenerate.
826    pub fn transformed(&self, t: &Transform2, tol: Tolerances) -> OgeomResult<Self> {
827        let scale = t.scale_factor().abs();
828        Ok(match self {
829            Self::Line(c) => Self::Line(Line2d {
830                axis: Axis2::new(
831                    t.apply(c.axis.location),
832                    Direction2::new(t.apply_vector(c.axis.direction.vector()), tol)?,
833                ),
834                // A line's parameter is a length, so the domain rescales with it.
835                domain: (c.domain.0 * scale, c.domain.1 * scale),
836            }),
837            Self::Circle(c) => Self::Circle(Circle2d {
838                circle: c.circle.transformed(t, tol)?,
839                ..*c
840            }),
841            Self::Ellipse(c) => Self::Ellipse(Ellipse2d {
842                ellipse: c.ellipse.transformed(t, tol)?,
843                ..*c
844            }),
845            Self::BSpline(c) => Self::BSpline(BSpline2d {
846                control: c
847                    .control
848                    .iter()
849                    .map(|w| Weighted::new(t.apply(w.point()), w.weight, tol))
850                    .collect::<OgeomResult<Vec<_>>>()?,
851                ..c.clone()
852            }),
853            Self::Offset(c) => Self::Offset(Box::new(Offset2d {
854                basis: c.basis.transformed(t, tol)?,
855                distance: c.distance * scale,
856            })),
857            Self::Trig(c) => Self::Trig(Trig2d {
858                c: t.apply(c.c),
859                d: t.apply_vector(c.d),
860                a: t.apply_vector(c.a),
861                b: t.apply_vector(c.b),
862                ..*c
863            }),
864            Self::Trimmed(c) => Self::Trimmed(Box::new(Trimmed2d {
865                basis: c.basis.transformed(t, tol)?,
866                domain: if matches!(c.basis, Self::Line(_)) {
867                    (c.domain.0 * scale, c.domain.1 * scale)
868                } else {
869                    c.domain
870                },
871                reversed: c.reversed,
872            })),
873        })
874    }
875}
876
877impl Reversible for PlanarCurve {
878    fn reversed(&self) -> Self {
879        match self {
880            Self::Line(c) => Self::Line(Line2d {
881                axis: Axis2::new(
882                    c.axis.point_at(c.domain.0 + c.domain.1),
883                    c.axis.direction.reversed(),
884                ),
885                domain: c.domain,
886            }),
887            Self::Circle(c) => Self::Circle(Circle2d {
888                reversed: !c.reversed,
889                ..*c
890            }),
891            Self::Ellipse(c) => Self::Ellipse(Ellipse2d {
892                reversed: !c.reversed,
893                ..*c
894            }),
895            Self::BSpline(c) => {
896                let (knots, control) = bspline::reverse(&c.knots, &c.control);
897                Self::BSpline(BSpline2d {
898                    knots,
899                    control,
900                    ..c.clone()
901                })
902            }
903            Self::Trimmed(c) => Self::Trimmed(Box::new(Trimmed2d {
904                reversed: !c.reversed,
905                ..(**c).clone()
906            })),
907            // Reversing flips the tangent and with it the left normal, so
908            // the distance negates to keep the same point set traversed
909            // backwards.
910            Self::Offset(c) => Self::Offset(Box::new(Offset2d {
911                basis: c.basis.reversed(),
912                distance: -c.distance,
913            })),
914            Self::Trig(c) => Self::Trig(Trig2d {
915                reversed: !c.reversed,
916                ..*c
917            }),
918        }
919    }
920}
921
922/// The angle a planar curve's tangent makes with the `u` axis at `t`.
923///
924/// The natural way to ask which way a pcurve is heading, which is what wire
925/// ordering and outer/inner classification are built from.
926///
927/// # Errors
928///
929/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) at a cusp, where
930/// the derivative vanishes.
931pub fn tangent_angle(curve: &PlanarCurve, t: f64, tol: Tolerances) -> OgeomResult<f64> {
932    let d = curve.d1_at(t, tol)?;
933    if d.is_zero(tol) {
934        ogeom_bail!(Construction, "curve has no tangent at {t}");
935    }
936    Ok(elementary::wrap_signed_angle(d.y.atan2(d.x)))
937}
938
939impl From<Line2d> for PlanarCurve {
940    fn from(c: Line2d) -> Self {
941        Self::Line(c)
942    }
943}
944impl From<Circle2d> for PlanarCurve {
945    fn from(c: Circle2d) -> Self {
946        Self::Circle(c)
947    }
948}
949impl From<Ellipse2d> for PlanarCurve {
950    fn from(c: Ellipse2d) -> Self {
951        Self::Ellipse(c)
952    }
953}
954impl From<BSpline2d> for PlanarCurve {
955    fn from(c: BSpline2d) -> Self {
956        Self::BSpline(c)
957    }
958}
959impl From<Trimmed2d> for PlanarCurve {
960    fn from(c: Trimmed2d) -> Self {
961        Self::Trimmed(Box::new(c))
962    }
963}
964
965#[cfg(test)]
966#[allow(clippy::unwrap_used)]
967mod tests {
968    use super::*;
969    use approx::assert_relative_eq;
970    use ogeom_math::Frame2;
971
972    #[test]
973    fn an_offset_2d_circle_is_the_larger_circle() {
974        use ogeom_math::Circle2;
975        // Counterclockwise circle: the left normal points inward, so a
976        // negative distance grows the radius.
977        let frame = Frame2::new(Point2::new(1.0, -2.0), Direction2::X);
978        let circle = PlanarCurve::Circle(Circle2d::new(Circle2::new(frame, 2.0, T).unwrap()));
979        let offset = Offset2d::new(circle, -1.0).unwrap();
980        for i in 0..8 {
981            let t = core::f64::consts::TAU * f64::from(i) / 8.0;
982            let p = offset.point_at(t, T).unwrap();
983            assert!((p.distance(Point2::new(1.0, -2.0)) - 3.0).abs() < 1e-12);
984        }
985        let h = 1e-6;
986        let d = offset.d1_at(1.0, T).unwrap();
987        let fd = (offset.point_at(1.0 + h, T).unwrap() - offset.point_at(1.0 - h, T).unwrap())
988            / (2.0 * h);
989        assert!((d - fd).magnitude() < 1e-5);
990        assert!(offset.derivatives_at(1.0, 2, T).is_err());
991    }
992
993    const T: Tolerances = Tolerances::millimetres();
994
995    fn frame() -> Frame2 {
996        Frame2::new(Point2::new(2.0, -1.0), Direction2::from_angle(0.4))
997    }
998
999    fn every_curve() -> Vec<PlanarCurve> {
1000        let spline = {
1001            let control = vec![
1002                Point2::new(0.0, 0.0),
1003                Point2::new(1.0, 2.0),
1004                Point2::new(3.0, 1.0),
1005                Point2::new(5.0, 0.0),
1006                Point2::new(6.0, -1.0),
1007            ];
1008            BSpline2d::new(
1009                KnotVector::clamped_uniform(3, control.len()).unwrap(),
1010                control,
1011                T,
1012            )
1013            .unwrap()
1014        };
1015        vec![
1016            Line2d::segment(Point2::ORIGIN, Point2::new(3.0, 4.0), T)
1017                .unwrap()
1018                .into(),
1019            Circle2d::new(Circle2::new(frame(), 2.0, T).unwrap()).into(),
1020            Ellipse2d::new(Ellipse2::new(frame(), 5.0, 3.0, T).unwrap()).into(),
1021            spline.clone().into(),
1022            Trimmed2d::new(spline.into(), 0.2, 0.8, T).unwrap().into(),
1023        ]
1024    }
1025
1026    fn interior(c: &PlanarCurve, n: usize) -> Vec<f64> {
1027        let (a, b) = c.domain();
1028        (1..n)
1029            .map(|i| {
1030                #[allow(clippy::cast_precision_loss)]
1031                let t = i as f64 / n as f64 + 0.0413;
1032                a + (b - a) * t
1033            })
1034            .collect()
1035    }
1036
1037    #[test]
1038    fn derivatives_agree_with_finite_differences() {
1039        let h = 1e-6;
1040        for c in every_curve() {
1041            for u in interior(&c, 8) {
1042                let d1 = c.d1_at(u, T).unwrap();
1043                let numeric = (c.point_at(u + h, T).unwrap() - c.point_at(u - h, T).unwrap())
1044                    * (1.0 / (2.0 * h));
1045                assert!(
1046                    (d1 - numeric).magnitude() <= 1e-5 * numeric.magnitude().max(1.0),
1047                    "{:?} at {u}",
1048                    c.kind()
1049                );
1050            }
1051        }
1052    }
1053
1054    #[test]
1055    fn derivatives_at_zero_returns_the_point() {
1056        for c in every_curve() {
1057            for u in interior(&c, 4) {
1058                let d = c.derivatives_at(u, 0, T).unwrap();
1059                assert_eq!(d.len(), 1);
1060                assert!(Point2::from_vector(d[0]).is_equal(c.point_at(u, T).unwrap(), T));
1061            }
1062        }
1063    }
1064
1065    #[test]
1066    fn out_of_domain_parameters_follow_periodicity() {
1067        for c in every_curve() {
1068            let (a, b) = c.domain();
1069            if c.is_periodic() {
1070                assert!(c.point_at(b + 1.0, T).is_ok(), "{:?}", c.kind());
1071            } else {
1072                assert!(c.point_at(b + 1.0, T).is_err(), "{:?}", c.kind());
1073                assert!(c.point_at(a - 1.0, T).is_err(), "{:?}", c.kind());
1074            }
1075        }
1076    }
1077
1078    #[test]
1079    fn reversal_traverses_the_same_points_backwards() {
1080        for c in every_curve() {
1081            let r = c.reversed();
1082            let (a, b) = c.domain();
1083            assert_eq!(r.domain(), (a, b), "{:?} changed its domain", c.kind());
1084            for i in 0..=8 {
1085                let t = f64::from(i) / 8.0;
1086                let forward = c.point_at(a + (b - a) * t, T).unwrap();
1087                let backward = r.point_at(a + (b - a) * (1.0 - t), T).unwrap();
1088                assert!(forward.is_equal(backward, T), "{:?} at {t}", c.kind());
1089            }
1090        }
1091    }
1092
1093    #[test]
1094    fn reversing_twice_is_the_identity() {
1095        for c in every_curve() {
1096            let twice = c.reversed().reversed();
1097            for u in interior(&c, 8) {
1098                assert!(
1099                    c.point_at(u, T)
1100                        .unwrap()
1101                        .is_equal(twice.point_at(u, T).unwrap(), T),
1102                    "{:?}",
1103                    c.kind()
1104                );
1105            }
1106        }
1107    }
1108
1109    #[test]
1110    fn a_reversed_curve_heads_the_other_way() {
1111        for c in every_curve() {
1112            let r = c.reversed();
1113            let (a, b) = c.domain();
1114            let u = a + (b - a) * 0.4;
1115            let forward = tangent_angle(&c, u, T).unwrap();
1116            let backward = tangent_angle(&r, mirror(u, a, b), T).unwrap();
1117            let difference = (forward - backward).abs();
1118            assert!(
1119                (difference - core::f64::consts::PI).abs() < 1e-9,
1120                "{:?}: {forward} vs {backward}",
1121                c.kind()
1122            );
1123        }
1124    }
1125
1126    #[test]
1127    fn transforms_move_curves() {
1128        let t = Transform2::rotation(Point2::new(1.0, 1.0), 0.7);
1129        for c in every_curve() {
1130            let moved = c.transformed(&t, T).unwrap();
1131            assert_eq!(moved.kind(), c.kind());
1132            for u in interior(&c, 8) {
1133                let expected = t.apply(c.point_at(u, T).unwrap());
1134                assert!(
1135                    moved.point_at(u, T).unwrap().is_equal(expected, T),
1136                    "{:?} at {u}",
1137                    c.kind()
1138                );
1139            }
1140        }
1141    }
1142
1143    #[test]
1144    fn a_line_segments_parameter_is_arc_length() {
1145        let l = Line2d::segment(Point2::ORIGIN, Point2::new(3.0, 4.0), T).unwrap();
1146        assert_eq!(l.domain(), (0.0, 5.0));
1147        assert!(
1148            l.point_at(2.5, T)
1149                .unwrap()
1150                .is_equal(Point2::new(1.5, 2.0), T)
1151        );
1152        assert_relative_eq!(l.d1_at(1.0, T).unwrap().magnitude(), 1.0, epsilon = 1e-15);
1153        assert!(Line2d::segment(Point2::ORIGIN, Point2::ORIGIN, T).is_err());
1154        assert!(Line2d::over(Axis2::X, 1.0, 1.0).is_err());
1155    }
1156
1157    #[test]
1158    fn a_rational_quadratic_traces_an_exact_arc() {
1159        let w = core::f64::consts::FRAC_1_SQRT_2;
1160        let control: Vec<_> = [
1161            (Point2::new(1.0, 0.0), 1.0),
1162            (Point2::new(1.0, 1.0), w),
1163            (Point2::new(0.0, 1.0), 1.0),
1164        ]
1165        .iter()
1166        .map(|(p, w)| Weighted::new(*p, *w, T).unwrap())
1167        .collect();
1168        let c = BSpline2d::rational(KnotVector::clamped_uniform(2, 3).unwrap(), control).unwrap();
1169        assert!(c.is_rational());
1170        for i in 0..=20 {
1171            let u = f64::from(i) / 20.0;
1172            assert_relative_eq!(
1173                c.point_at(u, T).unwrap().to_vector().magnitude(),
1174                1.0,
1175                epsilon = 1e-14
1176            );
1177        }
1178    }
1179
1180    #[test]
1181    fn tangent_angle_reports_the_heading() {
1182        // Along +x, then +y after a quarter turn of the circle.
1183        let l: PlanarCurve = Line2d::segment(Point2::ORIGIN, Point2::new(5.0, 0.0), T)
1184            .unwrap()
1185            .into();
1186        assert_relative_eq!(tangent_angle(&l, 1.0, T).unwrap(), 0.0, epsilon = 1e-15);
1187
1188        let c: PlanarCurve =
1189            Circle2d::new(Circle2::centred(Point2::ORIGIN, 1.0, T).unwrap()).into();
1190        assert_relative_eq!(
1191            tangent_angle(&c, 0.0, T).unwrap(),
1192            core::f64::consts::FRAC_PI_2,
1193            epsilon = 1e-12
1194        );
1195        assert_relative_eq!(
1196            tangent_angle(&c, core::f64::consts::FRAC_PI_2, T).unwrap(),
1197            core::f64::consts::PI,
1198            epsilon = 1e-12
1199        );
1200    }
1201
1202    #[test]
1203    fn trimming_is_bounds_checked_and_agrees_with_its_basis() {
1204        let base: PlanarCurve = Line2d::over(Axis2::X, 0.0, 10.0).unwrap().into();
1205        assert!(Trimmed2d::new(base.clone(), 2.0, 8.0, T).is_ok());
1206        assert!(Trimmed2d::new(base.clone(), 8.0, 2.0, T).is_err());
1207        assert!(Trimmed2d::new(base.clone(), -1.0, 5.0, T).is_err());
1208
1209        let trimmed = Trimmed2d::new(base.clone(), 2.0, 8.0, T).unwrap();
1210        assert_eq!(trimmed.domain(), (2.0, 8.0));
1211        for i in 0..=6 {
1212            let u = 2.0 + 6.0 * f64::from(i) / 6.0;
1213            assert!(
1214                trimmed
1215                    .point_at(u, T)
1216                    .unwrap()
1217                    .is_equal(base.point_at(u, T).unwrap(), T)
1218            );
1219        }
1220        assert!(trimmed.point_at(1.0, T).is_err());
1221    }
1222
1223    #[test]
1224    fn a_circle_in_parameter_space_closes_on_itself() {
1225        let c: PlanarCurve = Circle2d::new(Circle2::new(frame(), 2.0, T).unwrap()).into();
1226        assert!(c.is_closed(T) && c.is_periodic());
1227        let base = c.point_at(0.7, T).unwrap();
1228        for k in [-2.0_f64, 1.0, 3.0] {
1229            assert!(base.is_equal(c.point_at(k.mul_add(TAU, 0.7), T).unwrap(), T));
1230        }
1231    }
1232}
ogeom_geom/curve2d.rs

ogeom_geom/
curve2d.rs
