dotloom_geometry/
curve.rs

1//! Elementary curve pieces: line segments, circles, arcs and Bézier curves.
2
3use core::f64::consts::{PI, TAU};
4
5use serde::{Deserialize, Serialize};
6
7use crate::{
8    Aabb, Affine, GeoResult, GeometryError, LinearKind, Orientation, Point, Vector, error::finite, normalize_angle,
9    orientation,
10};
11
12/// Relative tolerance used to classify transforms as similarities.
13pub(crate) const SIMILARITY_REL: f64 = 1e-9;
14/// Recursion cap for adaptive subdivision (2^16 pieces per curve at most).
15const MAX_DEPTH: u32 = 16;
16
17/// A straight line segment from `a` to `b`.
18#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
19pub struct Segment {
20    /// Start point.
21    pub a: Point,
22    /// End point.
23    pub b: Point,
24}
25
26impl Segment {
27    /// Create a segment.
28    #[must_use]
29    pub const fn new(a: Point, b: Point) -> Self {
30        Self { a, b }
31    }
32
33    /// Vector `b - a`.
34    #[must_use]
35    pub fn vector(self) -> Vector {
36        self.b - self.a
37    }
38
39    /// Length.
40    #[must_use]
41    pub fn length(self) -> f64 {
42        self.a.distance(self.b)
43    }
44
45    /// Unit direction, `None` when degenerate.
46    #[must_use]
47    pub fn direction(self) -> Option<Vector> {
48        self.vector().normalize()
49    }
50
51    /// Point at parameter `t` (`0 → a`, `1 → b`).
52    #[must_use]
53    pub fn point_at(self, t: f64) -> Point {
54        self.a.lerp(self.b, t)
55    }
56
57    /// Midpoint.
58    #[must_use]
59    pub fn midpoint(self) -> Point {
60        self.a.midpoint(self.b)
61    }
62
63    /// Unclamped parameter of the orthogonal projection of `p` on the line.
64    /// Degenerate segments return `0`.
65    #[must_use]
66    pub fn line_param(self, p: Point) -> f64 {
67        let v = self.vector();
68        let len2 = v.length_sq();
69        if len2 > 0.0 && len2.is_finite() { (p - self.a).dot(v) / len2 } else { 0.0 }
70    }
71
72    /// Closest point on the segment and its parameter.
73    #[must_use]
74    pub fn closest(self, p: Point) -> (f64, Point) {
75        let t = self.line_param(p).clamp(0.0, 1.0);
76        (t, self.point_at(t))
77    }
78
79    /// Distance from `p` to the segment.
80    #[must_use]
81    pub fn distance_to_point(self, p: Point) -> f64 {
82        self.closest(p).1.distance(p)
83    }
84
85    /// Bounding box.
86    #[must_use]
87    pub fn bbox(self) -> Aabb {
88        Aabb::from_corners(self.a, self.b)
89    }
90
91    /// Reversed segment.
92    #[must_use]
93    pub const fn reversed(self) -> Self {
94        Self::new(self.b, self.a)
95    }
96
97    /// Transformed segment (any affine transform is exact for segments).
98    #[must_use]
99    pub fn transform(self, t: Affine) -> Self {
100        Self::new(t.apply(self.a), t.apply(self.b))
101    }
102}
103
104/// A full circle.
105#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
106pub struct Circle {
107    /// Center.
108    pub center: Point,
109    /// Radius (> 0 for a valid circle).
110    pub radius: f64,
111}
112
113impl Circle {
114    /// Create a validated circle.
115    pub fn new(center: Point, radius: f64) -> GeoResult<Self> {
116        if !center.is_finite() {
117            return Err(GeometryError::NonFinite("circle center"));
118        }
119        finite(radius, "circle radius")?;
120        if radius <= 0.0 {
121            return Err(GeometryError::Degenerate("circle radius must be > 0"));
122        }
123        Ok(Self { center, radius })
124    }
125
126    /// Point at angle `a` (radians).
127    #[must_use]
128    pub fn point_at_angle(self, a: f64) -> Point {
129        self.center + Vector::from_angle(a) * self.radius
130    }
131
132    /// Circumference.
133    #[must_use]
134    pub fn length(self) -> f64 {
135        TAU * self.radius
136    }
137
138    /// Bounding box.
139    #[must_use]
140    pub fn bbox(self) -> Aabb {
141        Aabb::from_corners(self.center, self.center).inflate(self.radius.abs())
142    }
143
144    /// Closest point on the circle; the center maps to angle 0.
145    #[must_use]
146    pub fn closest(self, p: Point) -> (f64, Point) {
147        let a = match (p - self.center).normalize() {
148            Some(v) => normalize_angle(v.angle()),
149            None => 0.0,
150        };
151        (a, self.point_at_angle(a))
152    }
153
154    /// Distance from `p` to the circle line.
155    #[must_use]
156    pub fn distance_to_point(self, p: Point) -> f64 {
157        (p.distance(self.center) - self.radius).abs()
158    }
159
160    /// Apply a transform; only similarities keep a circle a circle.
161    pub fn transform(self, t: Affine) -> GeoResult<Self> {
162        match t.linear_kind(SIMILARITY_REL) {
163            LinearKind::Similarity { scale, .. } => Self::new(t.apply(self.center), self.radius * scale),
164            LinearKind::General => Err(GeometryError::UnsupportedTransform {
165                shape: "circle",
166                reason: "non-uniform scale or shear turns a circle into an ellipse",
167            }),
168            LinearKind::Singular => Err(GeometryError::SingularTransform { determinant: t.determinant() }),
169        }
170    }
171}
172
173/// A circular arc: `center + radius·(cos θ, sin θ)` for `θ = start + sweep·t`, `t ∈ [0, 1]`.
174///
175/// `sweep > 0` is counter-clockwise. `|sweep| ≤ 2π`.
176#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
177pub struct Arc {
178    /// Center.
179    pub center: Point,
180    /// Radius (> 0).
181    pub radius: f64,
182    /// Start angle in radians.
183    pub start: f64,
184    /// Signed sweep in radians.
185    pub sweep: f64,
186}
187
188impl Arc {
189    /// Create a validated arc.
190    pub fn new(center: Point, radius: f64, start: f64, sweep: f64) -> GeoResult<Self> {
191        if !center.is_finite() {
192            return Err(GeometryError::NonFinite("arc center"));
193        }
194        finite(radius, "arc radius")?;
195        finite(start, "arc start angle")?;
196        finite(sweep, "arc sweep")?;
197        if radius <= 0.0 {
198            return Err(GeometryError::Degenerate("arc radius must be > 0"));
199        }
200        if sweep == 0.0 {
201            return Err(GeometryError::Degenerate("arc sweep must be non-zero"));
202        }
203        Ok(Self { center, radius, start, sweep: sweep.clamp(-TAU, TAU) })
204    }
205
206    /// Arc through three points `a → m → b`.
207    pub fn from_three_points(a: Point, m: Point, b: Point) -> GeoResult<Self> {
208        if !(a.is_finite() && m.is_finite() && b.is_finite()) {
209            return Err(GeometryError::NonFinite("arc points"));
210        }
211        let o = orientation(a, m, b);
212        if o == Orientation::Collinear {
213            return Err(GeometryError::Degenerate("arc points are collinear"));
214        }
215        let center = circumcenter(a, m, b).ok_or(GeometryError::Degenerate("arc points are collinear"))?;
216        let radius = center.distance(a);
217        let sa = (a - center).angle();
218        let sb = (b - center).angle();
219        let ccw = o == Orientation::CounterClockwise;
220        let mut sweep = normalize_angle(sb - sa);
221        if !ccw {
222            sweep -= TAU;
223        }
224        if sweep == 0.0 {
225            return Err(GeometryError::Degenerate("arc start and end coincide"));
226        }
227        Self::new(center, radius, sa, sweep)
228    }
229
230    /// Arc between `p0` and `p1` with DXF-style bulge (`tan(sweep/4)`, positive = CCW).
231    pub fn from_bulge(p0: Point, p1: Point, bulge: f64) -> GeoResult<Self> {
232        finite(bulge, "bulge")?;
233        if bulge == 0.0 {
234            return Err(GeometryError::Degenerate("zero bulge is a straight segment"));
235        }
236        let chord = p1 - p0;
237        let c = chord.length();
238        if c == 0.0 || !c.is_finite() {
239            return Err(GeometryError::Degenerate("bulge arc with coincident endpoints"));
240        }
241        let sweep = 4.0 * crate::math::atan(bulge);
242        let left = chord.perp() / c;
243        let h = (c * 0.5) * (1.0 - bulge * bulge) / (2.0 * bulge);
244        let center = p0.midpoint(p1) + left * h;
245        let radius = center.distance(p0);
246        Self::new(center, radius, (p0 - center).angle(), sweep)
247    }
248
249    /// Bulge value of this arc (`tan(sweep/4)`).
250    #[must_use]
251    pub fn bulge(self) -> f64 {
252        crate::math::tan(self.sweep / 4.0)
253    }
254
255    /// End angle (`start + sweep`).
256    #[must_use]
257    pub fn end_angle(self) -> f64 {
258        self.start + self.sweep
259    }
260
261    /// Point at parameter `t`.
262    #[must_use]
263    pub fn point_at(self, t: f64) -> Point {
264        self.center + Vector::from_angle(self.start + self.sweep * t) * self.radius
265    }
266
267    /// Start point.
268    #[must_use]
269    pub fn start_point(self) -> Point {
270        self.point_at(0.0)
271    }
272
273    /// End point.
274    #[must_use]
275    pub fn end_point(self) -> Point {
276        self.point_at(1.0)
277    }
278
279    /// Midpoint along the arc.
280    #[must_use]
281    pub fn mid_point(self) -> Point {
282        self.point_at(0.5)
283    }
284
285    /// Unit tangent at `t` in the direction of travel.
286    #[must_use]
287    pub fn tangent_at(self, t: f64) -> Vector {
288        let v = Vector::from_angle(self.start + self.sweep * t).perp();
289        if self.sweep >= 0.0 { v } else { -v }
290    }
291
292    /// Arc length.
293    #[must_use]
294    pub fn length(self) -> f64 {
295        self.radius * self.sweep.abs()
296    }
297
298    /// Parameter for absolute angle `a` if it lies on the arc (with `eps` radians slack).
299    #[must_use]
300    pub fn param_of_angle(self, a: f64, eps: f64) -> Option<f64> {
301        let d = if self.sweep >= 0.0 { normalize_angle(a - self.start) } else { normalize_angle(self.start - a) };
302        let span = self.sweep.abs();
303        if d <= span + eps {
304            Some((d / span).min(1.0))
305        } else if TAU - d <= eps {
306            Some(0.0)
307        } else {
308            None
309        }
310    }
311
312    /// Closest point on the arc and its parameter.
313    #[must_use]
314    pub fn closest(self, p: Point) -> (f64, Point) {
315        if let Some(v) = (p - self.center).normalize()
316            && let Some(t) = self.param_of_angle(v.angle(), 0.0)
317        {
318            return (t, self.point_at(t));
319        }
320        let (s, e) = (self.start_point(), self.end_point());
321        if p.distance_sq(s) <= p.distance_sq(e) { (0.0, s) } else { (1.0, e) }
322    }
323
324    /// Distance from `p` to the arc.
325    #[must_use]
326    pub fn distance_to_point(self, p: Point) -> f64 {
327        self.closest(p).1.distance(p)
328    }
329
330    /// Exact bounding box.
331    #[must_use]
332    pub fn bbox(self) -> Aabb {
333        let mut b = Aabb::from_corners(self.start_point(), self.end_point());
334        for k in 0..4 {
335            let a = f64::from(k) * PI * 0.5;
336            if self.param_of_angle(a, 0.0).is_some() {
337                b = b.include(self.center + Vector::from_angle(a) * self.radius);
338            }
339        }
340        b
341    }
342
343    /// Reversed arc (same points, opposite direction).
344    #[must_use]
345    pub fn reversed(self) -> Self {
346        Self { start: self.start + self.sweep, sweep: -self.sweep, ..self }
347    }
348
349    /// Sub-arc between parameters `t0` and `t1`.
350    #[must_use]
351    pub fn subarc(self, t0: f64, t1: f64) -> Self {
352        Self { start: self.start + self.sweep * t0, sweep: self.sweep * (t1 - t0), ..self }
353    }
354
355    /// Apply a similarity transform.
356    pub fn transform(self, t: Affine) -> GeoResult<Self> {
357        match t.linear_kind(SIMILARITY_REL) {
358            LinearKind::Similarity { scale, reflected } => {
359                let start_dir = t.apply_vector(Vector::from_angle(self.start));
360                let sweep = if reflected { -self.sweep } else { self.sweep };
361                Self::new(t.apply(self.center), self.radius * scale, start_dir.angle(), sweep)
362            }
363            LinearKind::General => Err(GeometryError::UnsupportedTransform {
364                shape: "arc",
365                reason: "non-uniform scale or shear turns an arc into an elliptical arc",
366            }),
367            LinearKind::Singular => Err(GeometryError::SingularTransform { determinant: t.determinant() }),
368        }
369    }
370
371    /// Append a polyline approximation (excluding the start point) to `out`.
372    pub fn flatten_into(self, tol: f64, out: &mut Vec<Point>) {
373        let n = arc_segments(self.radius, self.sweep.abs(), tol);
374        for i in 1..=n {
375            out.push(self.point_at(f64::from(i) / f64::from(n)));
376        }
377    }
378}
379
380/// Number of chords so that the sagitta stays below `tol`.
381pub(crate) fn arc_segments(radius: f64, sweep: f64, tol: f64) -> u32 {
382    let tol = tol.max(radius * 1e-12).max(1e-12);
383    if radius <= tol {
384        return 1;
385    }
386    // sagitta = r (1 - cos(θ/2)) ≤ tol  ⇒  θ ≤ 2 acos(1 - tol/r)
387    let max_step = 2.0 * crate::math::acos((1.0 - tol / radius).clamp(-1.0, 1.0));
388    if max_step <= 0.0 || !max_step.is_finite() {
389        return 4096;
390    }
391    let n = (sweep / max_step).ceil();
392    if n.is_finite() { n.clamp(1.0, 4096.0) as u32 } else { 4096 }
393}
394
395/// Circumcenter of a triangle, `None` when collinear.
396#[must_use]
397pub fn circumcenter(a: Point, b: Point, c: Point) -> Option<Point> {
398    let ab = b - a;
399    let ac = c - a;
400    let d = 2.0 * ab.cross(ac);
401    if d == 0.0 || !d.is_finite() {
402        return None;
403    }
404    let ab2 = ab.length_sq();
405    let ac2 = ac.length_sq();
406    let ux = (ac.y * ab2 - ab.y * ac2) / d;
407    let uy = (ab.x * ac2 - ac.x * ab2) / d;
408    let p = a + Vector::new(ux, uy);
409    p.is_finite().then_some(p)
410}
411
412/// Quadratic Bézier curve.
413#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
414pub struct QuadBez {
415    /// Start point.
416    pub p0: Point,
417    /// Control point.
418    pub p1: Point,
419    /// End point.
420    pub p2: Point,
421}
422
423/// Cubic Bézier curve.
424#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
425pub struct CubicBez {
426    /// Start point.
427    pub p0: Point,
428    /// First control point.
429    pub p1: Point,
430    /// Second control point.
431    pub p2: Point,
432    /// End point.
433    pub p3: Point,
434}
435
436impl QuadBez {
437    /// Elevate to an exactly equivalent cubic.
438    #[must_use]
439    pub fn to_cubic(self) -> CubicBez {
440        let c1 = self.p0 + (self.p1 - self.p0) * (2.0 / 3.0);
441        let c2 = self.p2 + (self.p1 - self.p2) * (2.0 / 3.0);
442        CubicBez { p0: self.p0, p1: c1, p2: c2, p3: self.p2 }
443    }
444}
445
446impl CubicBez {
447    /// Evaluate at `t`.
448    #[must_use]
449    pub fn eval(self, t: f64) -> Point {
450        let mt = 1.0 - t;
451        let a = mt * mt * mt;
452        let b = 3.0 * mt * mt * t;
453        let c = 3.0 * mt * t * t;
454        let d = t * t * t;
455        Point::new(
456            a * self.p0.x + b * self.p1.x + c * self.p2.x + d * self.p3.x,
457            a * self.p0.y + b * self.p1.y + c * self.p2.y + d * self.p3.y,
458        )
459    }
460
461    /// First derivative at `t`.
462    #[must_use]
463    pub fn deriv(self, t: f64) -> Vector {
464        let mt = 1.0 - t;
465        let d0 = self.p1 - self.p0;
466        let d1 = self.p2 - self.p1;
467        let d2 = self.p3 - self.p2;
468        d0 * (3.0 * mt * mt) + d1 * (6.0 * mt * t) + d2 * (3.0 * t * t)
469    }
470
471    /// Second derivative at `t`.
472    #[must_use]
473    pub fn deriv2(self, t: f64) -> Vector {
474        let dd0 = (self.p2 - self.p1) - (self.p1 - self.p0);
475        let dd1 = (self.p3 - self.p2) - (self.p2 - self.p1);
476        dd0 * (6.0 * (1.0 - t)) + dd1 * (6.0 * t)
477    }
478
479    /// Split at `t` (de Casteljau).
480    #[must_use]
481    pub fn split(self, t: f64) -> (Self, Self) {
482        let p01 = self.p0.lerp(self.p1, t);
483        let p12 = self.p1.lerp(self.p2, t);
484        let p23 = self.p2.lerp(self.p3, t);
485        let p012 = p01.lerp(p12, t);
486        let p123 = p12.lerp(p23, t);
487        let m = p012.lerp(p123, t);
488        (Self { p0: self.p0, p1: p01, p2: p012, p3: m }, Self { p0: m, p1: p123, p2: p23, p3: self.p3 })
489    }
490
491    /// Sub-curve between `t0` and `t1`.
492    #[must_use]
493    pub fn subsegment(self, t0: f64, t1: f64) -> Self {
494        if t1 <= t0 {
495            return Self { p0: self.eval(t0), p1: self.eval(t0), p2: self.eval(t0), p3: self.eval(t0) };
496        }
497        let (_, right) = self.split(t0);
498        let u = if t0 < 1.0 { (t1 - t0) / (1.0 - t0) } else { 1.0 };
499        right.split(u.clamp(0.0, 1.0)).0
500    }
501
502    /// Exact bounding box (extrema of the derivative).
503    #[must_use]
504    pub fn bbox(self) -> Aabb {
505        let mut b = Aabb::from_corners(self.p0, self.p3);
506        let roots_x = deriv_roots(self.p0.x, self.p1.x, self.p2.x, self.p3.x);
507        let roots_y = deriv_roots(self.p0.y, self.p1.y, self.p2.y, self.p3.y);
508        for t in roots_x.into_iter().chain(roots_y).flatten() {
509            if (0.0..=1.0).contains(&t) {
510                b = b.include(self.eval(t));
511            }
512        }
513        b
514    }
515
516    /// Bounding box of the control polygon (cheap, conservative).
517    #[must_use]
518    pub fn hull_bbox(self) -> Aabb {
519        Aabb::from_points([self.p0, self.p1, self.p2, self.p3])
520    }
521
522    /// Maximum distance of the control points from the chord (flatness bound).
523    #[must_use]
524    pub fn flatness(self) -> f64 {
525        let chord = Segment::new(self.p0, self.p3);
526        chord.distance_to_point(self.p1).max(chord.distance_to_point(self.p2))
527    }
528
529    /// Append a polyline approximation (excluding the start point) to `out`.
530    pub fn flatten_into(self, tol: f64, out: &mut Vec<Point>) {
531        let tol = tol.max(1e-12);
532        flatten_rec(self, tol, 0, out);
533    }
534
535    /// Arc length with adaptive Gauss–Legendre quadrature (relative accuracy ≈ `rel`).
536    #[must_use]
537    pub fn length(self, rel: f64) -> f64 {
538        length_rec(self, gauss_len(self), rel.max(1e-14), 0)
539    }
540
541    /// Closest point (parameter, point).
542    ///
543    /// Global search by branch-and-bound subdivision (the control-polygon box is a
544    /// lower bound for every sub-curve), followed by Newton polishing. This cannot
545    /// get trapped in a local minimum the way pure sampling + Newton can.
546    #[must_use]
547    pub fn closest(self, p: Point) -> (f64, Point) {
548        let size = self.hull_bbox().size().length();
549        let eps = (size * 1e-12).max(1e-300);
550        let mut best = (f64::INFINITY, 0.0);
551        for (t, q) in [(0.0, self.p0), (1.0, self.p3)] {
552            let d = q.distance_sq(p);
553            if d < best.0 {
554                best = (d, t);
555            }
556        }
557        closest_bb(self, 0.0, 1.0, p, eps, 0, &mut best);
558        let (best_d, best_t) = best;
559        let mut t = best_t;
560        for _ in 0..12 {
561            let b = self.eval(t);
562            let d1 = self.deriv(t);
563            let d2 = self.deriv2(t);
564            let f = (b - p).dot(d1);
565            let fp = d1.dot(d1) + (b - p).dot(d2);
566            if fp.abs() < 1e-300 || !fp.is_finite() {
567                break;
568            }
569            let nt = (t - f / fp).clamp(0.0, 1.0);
570            if (nt - t).abs() < 1e-15 {
571                t = nt;
572                break;
573            }
574            t = nt;
575        }
576        let refined = self.eval(t);
577        if refined.distance_sq(p) <= best_d { (t, refined) } else { (best_t, self.eval(best_t)) }
578    }
579
580    /// Transform (any affine transform is exact for Bézier curves).
581    #[must_use]
582    pub fn transform(self, t: Affine) -> Self {
583        Self { p0: t.apply(self.p0), p1: t.apply(self.p1), p2: t.apply(self.p2), p3: t.apply(self.p3) }
584    }
585
586    /// Reversed curve.
587    #[must_use]
588    pub const fn reversed(self) -> Self {
589        Self { p0: self.p3, p1: self.p2, p2: self.p1, p3: self.p0 }
590    }
591}
592
593fn closest_bb(c: CubicBez, t0: f64, t1: f64, p: Point, eps: f64, depth: u32, best: &mut (f64, f64)) {
594    let lb = c.hull_bbox().distance_to_point(p);
595    if lb * lb > best.0 {
596        return;
597    }
598    if depth >= 48 || c.flatness() <= eps {
599        let chord = Segment::new(c.p0, c.p3);
600        let (u, _) = chord.closest(p);
601        let t = t0 + (t1 - t0) * u;
602        let q = c.eval(u);
603        let d = q.distance_sq(p);
604        if d < best.0 {
605            *best = (d, t);
606        }
607        return;
608    }
609    let (l, r) = c.split(0.5);
610    let tm = 0.5 * (t0 + t1);
611    let mid = l.p3.distance_sq(p);
612    if mid < best.0 {
613        *best = (mid, tm);
614    }
615    // Visit the nearer half first for better pruning.
616    if l.hull_bbox().distance_to_point(p) <= r.hull_bbox().distance_to_point(p) {
617        closest_bb(l, t0, tm, p, eps, depth + 1, best);
618        closest_bb(r, tm, t1, p, eps, depth + 1, best);
619    } else {
620        closest_bb(r, tm, t1, p, eps, depth + 1, best);
621        closest_bb(l, t0, tm, p, eps, depth + 1, best);
622    }
623}
624
625fn flatten_rec(c: CubicBez, tol: f64, depth: u32, out: &mut Vec<Point>) {
626    if depth >= MAX_DEPTH || c.flatness() <= tol {
627        out.push(c.p3);
628        return;
629    }
630    let (l, r) = c.split(0.5);
631    flatten_rec(l, tol, depth + 1, out);
632    flatten_rec(r, tol, depth + 1, out);
633}
634
635/// Roots in [0,1] of the derivative of a 1D cubic Bézier.
636fn deriv_roots(p0: f64, p1: f64, p2: f64, p3: f64) -> [Option<f64>; 2] {
637    // B'(t)/3 = a t² + b t + c
638    let a = -p0 + 3.0 * p1 - 3.0 * p2 + p3;
639    let b = 2.0 * (p0 - 2.0 * p1 + p2);
640    let c = p1 - p0;
641    let scale = a.abs().max(b.abs()).max(c.abs());
642    if scale == 0.0 {
643        return [None, None];
644    }
645    if a.abs() <= 1e-12 * scale {
646        if b.abs() <= 1e-12 * scale {
647            return [None, None];
648        }
649        return [Some(-c / b), None];
650    }
651    let disc = b * b - 4.0 * a * c;
652    if disc < 0.0 {
653        return [None, None];
654    }
655    let sq = disc.sqrt();
656    // Numerically stable quadratic formula.
657    let q = -0.5 * (b + b.signum() * sq);
658    let r1 = q / a;
659    let r2 = if q != 0.0 { c / q } else { -b / (2.0 * a) };
660    [Some(r1), Some(r2)]
661}
662
663const GL_X: [f64; 8] = [
664    -0.960_289_856_497_536_2,
665    -0.796_666_477_413_626_7,
666    -0.525_532_409_916_329,
667    -0.183_434_642_495_649_8,
668    0.183_434_642_495_649_8,
669    0.525_532_409_916_329,
670    0.796_666_477_413_626_7,
671    0.960_289_856_497_536_2,
672];
673const GL_W: [f64; 8] = [
674    0.101_228_536_290_376_3,
675    0.222_381_034_453_374_5,
676    0.313_706_645_877_887_3,
677    0.362_683_783_378_362,
678    0.362_683_783_378_362,
679    0.313_706_645_877_887_3,
680    0.222_381_034_453_374_5,
681    0.101_228_536_290_376_3,
682];
683
684fn gauss_len(c: CubicBez) -> f64 {
685    GL_X.iter().zip(GL_W.iter()).map(|(x, w)| w * c.deriv(0.5 * (x + 1.0)).length()).sum::<f64>() * 0.5
686}
687
688fn length_rec(c: CubicBez, whole: f64, rel: f64, depth: u32) -> f64 {
689    let (l, r) = c.split(0.5);
690    let (ll, rl) = (gauss_len(l), gauss_len(r));
691    let halves = ll + rl;
692    if depth >= 12 || (halves - whole).abs() <= rel * halves.max(1e-300) {
693        halves
694    } else {
695        length_rec(l, ll, rel, depth + 1) + length_rec(r, rl, rel, depth + 1)
696    }
697}
698
699/// One elementary curve piece. Shapes decompose into these for intersections,
700/// trimming, hit-testing and tessellation.
701#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
702#[serde(tag = "kind", rename_all = "camelCase")]
703pub enum Curve {
704    /// Straight segment.
705    Line(Segment),
706    /// Circular arc.
707    Arc(Arc),
708    /// Full circle (closed, parameter = angle / 2π).
709    Circle(Circle),
710    /// Cubic Bézier (quadratics are elevated exactly).
711    Cubic(CubicBez),
712}
713
714impl Curve {
715    /// Point at normalized parameter `t ∈ [0,1]`.
716    #[must_use]
717    pub fn point_at(&self, t: f64) -> Point {
718        match *self {
719            Self::Line(s) => s.point_at(t),
720            Self::Arc(a) => a.point_at(t),
721            Self::Circle(c) => c.point_at_angle(t * TAU),
722            Self::Cubic(c) => c.eval(t),
723        }
724    }
725
726    /// Start point (circles start at angle 0).
727    #[must_use]
728    pub fn start(&self) -> Point {
729        self.point_at(0.0)
730    }
731
732    /// End point.
733    #[must_use]
734    pub fn end(&self) -> Point {
735        self.point_at(1.0)
736    }
737
738    /// Whether the curve is closed on itself.
739    #[must_use]
740    pub fn is_closed(&self) -> bool {
741        matches!(self, Self::Circle(_))
742    }
743
744    /// Bounding box.
745    #[must_use]
746    pub fn bbox(&self) -> Aabb {
747        match *self {
748            Self::Line(s) => s.bbox(),
749            Self::Arc(a) => a.bbox(),
750            Self::Circle(c) => c.bbox(),
751            Self::Cubic(c) => c.bbox(),
752        }
753    }
754
755    /// Closest point `(t, point)`.
756    #[must_use]
757    pub fn closest(&self, p: Point) -> (f64, Point) {
758        match *self {
759            Self::Line(s) => s.closest(p),
760            Self::Arc(a) => a.closest(p),
761            Self::Circle(c) => {
762                let (a, q) = c.closest(p);
763                (a / TAU, q)
764            }
765            Self::Cubic(c) => c.closest(p),
766        }
767    }
768
769    /// Distance to `p`.
770    #[must_use]
771    pub fn distance_to_point(&self, p: Point) -> f64 {
772        self.closest(p).1.distance(p)
773    }
774
775    /// Length.
776    #[must_use]
777    pub fn length(&self) -> f64 {
778        match *self {
779            Self::Line(s) => s.length(),
780            Self::Arc(a) => a.length(),
781            Self::Circle(c) => c.length(),
782            Self::Cubic(c) => c.length(1e-12),
783        }
784    }
785
786    /// Append flattened points (excluding the start point) to `out`.
787    pub fn flatten_into(&self, tol: f64, out: &mut Vec<Point>) {
788        match *self {
789            Self::Line(s) => out.push(s.b),
790            Self::Arc(a) => a.flatten_into(tol, out),
791            Self::Circle(c) => {
792                if let Ok(a) = Arc::new(c.center, c.radius, 0.0, TAU) {
793                    a.flatten_into(tol, out);
794                }
795            }
796            Self::Cubic(c) => c.flatten_into(tol, out),
797        }
798    }
799
800    /// Unit tangent at `t` (None where the derivative vanishes).
801    #[must_use]
802    pub fn tangent_at(&self, t: f64) -> Option<Vector> {
803        match *self {
804            Self::Line(s) => s.direction(),
805            Self::Arc(a) => Some(a.tangent_at(t)),
806            Self::Circle(_) => Some(Vector::from_angle(t * TAU).perp()),
807            Self::Cubic(c) => c.deriv(t).normalize(),
808        }
809    }
810
811    /// Kind name for diagnostics.
812    #[must_use]
813    pub const fn kind_name(&self) -> &'static str {
814        match self {
815            Self::Line(_) => "line",
816            Self::Arc(_) => "arc",
817            Self::Circle(_) => "circle",
818            Self::Cubic(_) => "cubic",
819        }
820    }
821}
822
823#[cfg(test)]
824mod tests {
825    use super::*;
826
827    #[test]
828    fn arc_from_three_points_ccw_and_cw() {
829        let a = Arc::from_three_points(Point::new(1.0, 0.0), Point::new(0.0, 1.0), Point::new(-1.0, 0.0)).unwrap();
830        assert!(a.center.distance(Point::ORIGIN) < 1e-12);
831        assert!((a.radius - 1.0).abs() < 1e-12);
832        assert!((a.sweep - PI).abs() < 1e-12);
833        let b = Arc::from_three_points(Point::new(1.0, 0.0), Point::new(0.0, -1.0), Point::new(-1.0, 0.0)).unwrap();
834        assert!((b.sweep + PI).abs() < 1e-12);
835        assert!(b.mid_point().distance(Point::new(0.0, -1.0)) < 1e-12);
836        assert!(Arc::from_three_points(Point::ORIGIN, Point::new(1.0, 1.0), Point::new(2.0, 2.0)).is_err());
837    }
838
839    #[test]
840    fn bulge_roundtrip() {
841        for bulge in [0.25, 1.0, -0.5, 2.0, -3.0] {
842            let p0 = Point::new(2.0, 1.0);
843            let p1 = Point::new(5.0, -3.0);
844            let a = Arc::from_bulge(p0, p1, bulge).unwrap();
845            assert!(a.start_point().distance(p0) < 1e-9, "bulge {bulge}");
846            assert!(a.end_point().distance(p1) < 1e-9, "bulge {bulge}");
847            assert!((a.bulge() - bulge).abs() < 1e-9);
848        }
849        // Semicircle: bulge 1, center at chord midpoint.
850        let s = Arc::from_bulge(Point::ORIGIN, Point::new(2.0, 0.0), 1.0).unwrap();
851        assert!(s.center.distance(Point::new(1.0, 0.0)) < 1e-12);
852        // CCW from (0,0) to (2,0) goes below the chord.
853        assert!(s.mid_point().y < 0.0);
854    }
855
856    #[test]
857    fn arc_bbox_includes_extrema() {
858        let a = Arc::new(Point::ORIGIN, 1.0, -0.1, 0.2 + PI / 2.0).unwrap();
859        let b = a.bbox();
860        assert!((b.max.x - 1.0).abs() < 1e-12);
861        assert!((b.max.y - 1.0).abs() < 1e-12);
862    }
863
864    #[test]
865    fn circle_rejects_nonuniform_scale() {
866        let c = Circle::new(Point::ORIGIN, 2.0).unwrap();
867        assert!(matches!(c.transform(Affine::scale(2.0, 1.0)), Err(GeometryError::UnsupportedTransform { .. })));
868        let r = c.transform(Affine::scale(3.0, 3.0)).unwrap();
869        assert!((r.radius - 6.0).abs() < 1e-12);
870    }
871
872    #[test]
873    fn arc_reflection_flips_sweep() {
874        let a = Arc::new(Point::ORIGIN, 1.0, 0.0, PI / 2.0).unwrap();
875        let m = Affine::mirror(Point::ORIGIN, Vector::new(1.0, 0.0)).unwrap();
876        let r = a.transform(m).unwrap();
877        assert!(r.end_point().distance(Point::new(0.0, -1.0)) < 1e-12);
878        assert!(r.sweep < 0.0);
879    }
880
881    #[test]
882    fn cubic_length_of_straight_line() {
883        let c = CubicBez {
884            p0: Point::ORIGIN,
885            p1: Point::new(1.0, 0.0),
886            p2: Point::new(2.0, 0.0),
887            p3: Point::new(3.0, 0.0),
888        };
889        assert!((c.length(1e-12) - 3.0).abs() < 1e-10);
890    }
891
892    #[test]
893    fn quarter_circle_cubic_length() {
894        // Standard approximation constant; true arc length differs by < 1e-3.
895        let k = 0.552_284_749_830_793_4;
896        let c = CubicBez {
897            p0: Point::new(1.0, 0.0),
898            p1: Point::new(1.0, k),
899            p2: Point::new(k, 1.0),
900            p3: Point::new(0.0, 1.0),
901        };
902        assert!((c.length(1e-12) - PI / 2.0).abs() < 1e-3);
903    }
904
905    #[test]
906    fn cubic_closest_matches_dense_sampling() {
907        let c = CubicBez {
908            p0: Point::ORIGIN,
909            p1: Point::new(1.0, 3.0),
910            p2: Point::new(4.0, -2.0),
911            p3: Point::new(5.0, 1.0),
912        };
913        let p = Point::new(2.3, 1.7);
914        let (_, q) = c.closest(p);
915        let mut best = f64::INFINITY;
916        for i in 0..=100_000 {
917            best = best.min(c.eval(f64::from(i) / 100_000.0).distance(p));
918        }
919        assert!(q.distance(p) <= best + 1e-9);
920    }
921
922    #[test]
923    fn flatten_respects_tolerance() {
924        let a = Arc::new(Point::ORIGIN, 100.0, 0.0, PI).unwrap();
925        let mut pts = vec![a.start_point()];
926        a.flatten_into(0.01, &mut pts);
927        for w in pts.windows(2) {
928            let mid = w[0].midpoint(w[1]);
929            assert!(100.0 - mid.distance(Point::ORIGIN) <= 0.01 + 1e-12);
930        }
931    }
932}
dotloom_geometry/curve.rs

dotloom_geometry/
curve.rs
