axiolid_nurbs/
exact_curve_intersection.rs

1//! Exact intersections of analytic curves with elementary surfaces and
2//! with each other (#119).
3//!
4//! # Method
5//!
6//! A line `o + t d` or a conic `o + x a cos(theta) + y b sin(theta)` is
7//! written as a polynomial vector over one parameter: `t` for a line, the
8//! half-angle `w = tan(theta / 2)` for a conic, using
9//! `cos = (1 - w^2) / (1 + w^2)` and `sin = 2 w / (1 + w^2)`. Substituting
10//! that into the implicit equation of the other operand and clearing the
11//! positive denominator gives an integer polynomial in one variable, whose
12//! real roots are isolated exactly by Sturm sequences (`axiolid-exact`).
13//!
14//! Every decision is an exact sign on the given `f64` input:
15//!
16//! - **containment**: the polynomial is identically zero, so the curve
17//!   lies in the other operand;
18//! - **count and order**: exact root isolation;
19//! - **tangency**: a root's multiplicity, decided by exact signs of
20//!   derivatives at the root;
21//! - **the conic point at `theta = pi`**, where `w` is infinite: the
22//!   polynomial loses degree, and the lost degree is that point's
23//!   multiplicity.
24//!
25//! # Which surface
26//!
27//! Implicit equations are derived from the same parametrisations
28//! `axiolid-evaluate` uses (`o + x r cos u + y r sin u + z v` for a
29//! cylinder, and so on), with frame axes taken exactly as given, not
30//! assumed orthonormal. Local coordinates come from Cramer's rule, so the
31//! equations stay polynomial in the input. A cone's slope is the `f64`
32//! value `tan(semi_angle)`, the value the evaluator also uses; the answer
33//! is exact for that cone. A cone point counts only on the nappe the
34//! evaluator covers (`radius + v * slope >= 0`).
35//!
36//! Only returned parameters are exact. Points are rounded, for output.
37
38use axiolid_core::{Frame3, Point3, Vec3};
39use axiolid_curve::{Curve2, Curve3};
40use axiolid_evaluate::evaluate3;
41use axiolid_exact::{Arith, Dyadic, IntPoly, RealRoot};
42use axiolid_guarantees::Sign;
43use axiolid_surface::Surface;
44
45/// Why an exact curve intersection was not computed.
46#[derive(Debug, Clone, PartialEq, Eq)]
47#[non_exhaustive]
48pub enum ExactCurveRefusal {
49    /// The curve is not a line, circle or ellipse. B-spline curves go
50    /// through the certified numeric path instead.
51    UnsupportedCurve,
52    /// The surface is a B-spline surface or an unknown kind.
53    UnsupportedSurface,
54    /// An input coordinate, radius or angle is NaN or infinite.
55    NonFinite,
56    /// A frame is singular, a direction is zero, or a radius is not
57    /// positive: there is no curve or surface to intersect.
58    Degenerate,
59    /// The curve lies in the cone's quadric but crosses its apex, so only
60    /// part of it is on the modelled nappe: the overlap is a ray, which
61    /// is neither a finite point set nor the whole curve.
62    PartialOverlap,
63}
64
65/// Where on the first curve an intersection lies, exactly.
66#[derive(Debug, Clone, PartialEq, Eq)]
67#[non_exhaustive]
68pub enum ExactCurveParameter {
69    /// A line's own parameter `t` (point `origin + t * direction`).
70    Line(RealRoot),
71    /// A conic's half-angle parameter `w = tan(theta / 2)`.
72    HalfAngle(RealRoot),
73    /// A conic's point at `theta = pi`, where `w` is infinite.
74    Antipode,
75    /// A section-family curve's own parameter, isolated by certified
76    /// subdivision (ADR 0077) and refined to the last bits.
77    Certified(Isolated),
78}
79
80/// A parameter isolated by certified subdivision, held as the exact double
81/// it was refined to (so parameters compare exactly, bit for bit).
82#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
83pub struct Isolated(u64);
84
85impl Isolated {
86    /// The parameter `value`.
87    #[must_use]
88    pub fn new(value: f64) -> Self {
89        Self(value.to_bits())
90    }
91
92    /// The parameter as a double.
93    #[must_use]
94    pub fn value(self) -> f64 {
95        f64::from_bits(self.0)
96    }
97}
98
99impl ExactCurveParameter {
100    /// The curve parameter as a double, for output: `t` for a line, the
101    /// angle `theta` in `[0, 2 pi)` for a conic.
102    #[must_use]
103    pub fn approx(&self) -> f64 {
104        match self {
105            Self::Line(t) => t.approx(),
106            Self::HalfAngle(w) => {
107                let theta = 2.0 * w.approx().atan();
108                if theta < 0.0 {
109                    theta + std::f64::consts::TAU
110                } else {
111                    theta
112                }
113            }
114            Self::Antipode => std::f64::consts::PI,
115            Self::Certified(t) => t.value(),
116        }
117    }
118}
119
120/// One intersection point.
121#[derive(Debug, Clone, PartialEq)]
122#[non_exhaustive]
123pub struct ExactCurveHit {
124    /// Its parameter on the first curve, exactly.
125    pub parameter: ExactCurveParameter,
126    /// Its multiplicity: 1 for a transverse crossing, 2 or more where the
127    /// curve touches the other operand.
128    pub multiplicity: usize,
129    /// The point, rounded from the curve at the parameter (output only).
130    /// Planar results have `z = 0`.
131    pub point: Point3,
132}
133
134impl ExactCurveHit {
135    /// Whether the curve touches rather than crosses here.
136    #[must_use]
137    pub fn is_tangent(&self) -> bool {
138        self.multiplicity >= 2
139    }
140}
141
142/// The exact intersection of a curve with another operand.
143#[derive(Debug, Clone, PartialEq)]
144#[non_exhaustive]
145pub enum ExactCurveIntersection {
146    /// The whole curve lies in the other operand.
147    Contained,
148    /// Finitely many points, ordered along the first curve: by `t` for a
149    /// line, by `theta` in `[0, 2 pi)` for a conic. Empty when they miss.
150    Points(Vec<ExactCurveHit>),
151}
152
153/// The exact intersection of an analytic curve with an elementary surface.
154///
155/// Curves: line, circle, ellipse. Surfaces: plane, cylinder, elliptical
156/// cylinder, cone, sphere, torus.
157///
158/// # Errors
159///
160/// [`ExactCurveRefusal`] names the unsupported or malformed operand.
161pub fn exact_curve_surface_intersection(
162    curve: &Curve3,
163    surface: &Surface,
164) -> Result<ExactCurveIntersection, ExactCurveRefusal> {
165    // A traced section: along its cells, with certified root isolation
166    // (ADR 0077). The ADR 0076 graphs need their span, so they take
167    // `section_curve_surface_intersection`.
168    if let Curve3::ImplicitSection(section) = curve {
169        let span = axiolid_core::Interval::new(0.0, section.curve.end());
170        return crate::implicit_ops::section_curve_surface_intersection(curve, span, surface);
171    }
172    // A B-spline curve on a B-spline surface: hull pruning over the Bezier
173    // pieces of both, then Newton (ADR 0077).
174    if let (Curve3::BSpline(c), Surface::BSpline(s)) = (curve, surface) {
175        let hits = crate::pair_trace::spline_curve_surface_hits(c, s)
176            .ok_or(ExactCurveRefusal::UnsupportedCurve)?;
177        return Ok(ExactCurveIntersection::Points(
178            hits.into_iter()
179                .map(|(t, point)| ExactCurveHit {
180                    parameter: ExactCurveParameter::Certified(Isolated::new(t)),
181                    multiplicity: 1,
182                    point,
183                })
184                .collect(),
185        ));
186    }
187    // A B-spline surface: traced (ADR 0077) for lines and conics.
188    if matches!(surface, Surface::BSpline(_)) {
189        return crate::implicit_ops::conic_spline_intersection(curve, surface);
190    }
191    let param = Param::of(curve)?;
192    let locus = surface_locus(surface, &param)?;
193    let result = solve(curve, &param, &locus);
194    if matches!(result, ExactCurveIntersection::Contained) {
195        if let Some(condition) = &locus.condition {
196            if !nonnegative_everywhere(condition) {
197                return Err(ExactCurveRefusal::PartialOverlap);
198            }
199        }
200    }
201    Ok(result)
202}
203
204/// Whether an integer polynomial is `>= 0` on the whole real line (and so
205/// also at a conic's antipode, its formal leading term).
206///
207/// A sign change needs a root of odd multiplicity; without one the sign
208/// away from roots is the leading coefficient's.
209fn nonnegative_everywhere(c: &DPoly) -> bool {
210    let int = c.to_int();
211    if int.is_zero() {
212        return true;
213    }
214    let odd = int
215        .real_roots()
216        .iter()
217        .any(|root| multiplicity(root, &int) % 2 == 1);
218    !odd && c.lead_sign() != Sign::Negative
219}
220
221/// The exact intersection of two analytic space curves.
222///
223/// Parameters are those of `first`; swap the arguments for `second`'s.
224///
225/// # Errors
226///
227/// [`ExactCurveRefusal`] names the unsupported or malformed operand.
228pub fn exact_curve_curve_intersection3(
229    first: &Curve3,
230    second: &Curve3,
231) -> Result<ExactCurveIntersection, ExactCurveRefusal> {
232    let param = Param::of(first)?;
233    let locus = curve_locus(second, &param)?;
234    Ok(solve(first, &param, &locus))
235}
236
237/// The exact intersection of two analytic plane curves.
238///
239/// Parameters are those of `first`; points carry `z = 0`.
240///
241/// # Errors
242///
243/// [`ExactCurveRefusal`] names the unsupported or malformed operand.
244pub fn exact_curve_curve_intersection2(
245    first: &Curve2,
246    second: &Curve2,
247) -> Result<ExactCurveIntersection, ExactCurveRefusal> {
248    exact_curve_curve_intersection3(&lift(first)?, &lift(second)?)
249}
250
251// --- polynomials over dyadics ----------------------------------------------
252
253/// A polynomial with dyadic coefficients, lowest degree first.
254#[derive(Debug, Clone)]
255struct DPoly(Vec<Dyadic>);
256
257impl DPoly {
258    fn constant(c: Dyadic) -> Self {
259        Self(vec![c])
260    }
261
262    fn add(&self, other: &Self) -> Self {
263        let n = self.0.len().max(other.0.len());
264        Self(
265            (0..n)
266                .map(|i| match (self.0.get(i), other.0.get(i)) {
267                    (Some(a), Some(b)) => a.add(b),
268                    (Some(a), None) | (None, Some(a)) => a.clone(),
269                    (None, None) => Dyadic::zero(),
270                })
271                .collect(),
272        )
273    }
274
275    fn scale(&self, c: &Dyadic) -> Self {
276        Self(self.0.iter().map(|a| a.mul(c)).collect())
277    }
278
279    fn sub(&self, other: &Self) -> Self {
280        self.add(&other.scale(&dy(-1.0)))
281    }
282
283    fn mul(&self, other: &Self) -> Self {
284        if self.0.is_empty() || other.0.is_empty() {
285            return Self(Vec::new());
286        }
287        let mut out = vec![Dyadic::zero(); self.0.len() + other.0.len() - 1];
288        for (i, a) in self.0.iter().enumerate() {
289            for (j, b) in other.0.iter().enumerate() {
290                out[i + j] = out[i + j].add(&a.mul(b));
291            }
292        }
293        Self(out)
294    }
295
296    fn square(&self) -> Self {
297        self.mul(self)
298    }
299
300    fn to_int(&self) -> IntPoly {
301        IntPoly::from_dyadic(&self.0)
302    }
303
304    /// Sign of the highest non-zero coefficient (zero for the zero
305    /// polynomial).
306    fn lead_sign(&self) -> Sign {
307        self.0
308            .iter()
309            .rev()
310            .filter_map(Arith::sign)
311            .find(|s| *s != Sign::Zero)
312            .unwrap_or(Sign::Zero)
313    }
314}
315
316fn dy(value: f64) -> Dyadic {
317    Dyadic::from_f64(value)
318}
319
320type V3 = [Dyadic; 3];
321
322fn v3(v: Vec3) -> V3 {
323    [dy(v.x), dy(v.y), dy(v.z)]
324}
325
326fn dot(a: &V3, b: &V3) -> Dyadic {
327    a[0].mul(&b[0]).add(&a[1].mul(&b[1])).add(&a[2].mul(&b[2]))
328}
329
330fn cross(a: &V3, b: &V3) -> V3 {
331    [
332        a[1].mul(&b[2]).sub(&a[2].mul(&b[1])),
333        a[2].mul(&b[0]).sub(&a[0].mul(&b[2])),
334        a[0].mul(&b[1]).sub(&a[1].mul(&b[0])),
335    ]
336}
337
338fn is_zero(v: &Dyadic) -> bool {
339    v.sign() == Some(Sign::Zero)
340}
341
342fn finite(values: &[f64]) -> Result<(), ExactCurveRefusal> {
343    if values.iter().all(|v| v.is_finite()) {
344        Ok(())
345    } else {
346        Err(ExactCurveRefusal::NonFinite)
347    }
348}
349
350fn positive(values: &[f64]) -> Result<(), ExactCurveRefusal> {
351    finite(values)?;
352    if values.iter().all(|v| *v > 0.0) {
353        Ok(())
354    } else {
355        Err(ExactCurveRefusal::Degenerate)
356    }
357}
358
359fn frame_finite(f: &Frame3) -> Result<(), ExactCurveRefusal> {
360    finite(&[
361        f.origin.x, f.origin.y, f.origin.z, f.x.x, f.x.y, f.x.z, f.y.x, f.y.y, f.y.z, f.z.x, f.z.y,
362        f.z.z,
363    ])
364}
365
366// --- the first curve as a polynomial vector --------------------------------
367
368/// `X(s) = n(s) / w(s)` with `w > 0` for every real `s`.
369struct Param {
370    n: [DPoly; 3],
371    w: DPoly,
372    /// Degree of `w`: 0 for a line, 2 for a conic (half-angle form).
373    w_degree: usize,
374}
375
376impl Param {
377    fn of(curve: &Curve3) -> Result<Self, ExactCurveRefusal> {
378        match curve {
379            Curve3::Line(l) => {
380                finite(&[
381                    l.origin.x,
382                    l.origin.y,
383                    l.origin.z,
384                    l.direction.x,
385                    l.direction.y,
386                    l.direction.z,
387                ])?;
388                let (o, d) = (v3(l.origin), v3(l.direction));
389                if d.iter().all(is_zero) {
390                    return Err(ExactCurveRefusal::Degenerate);
391                }
392                Ok(Self {
393                    n: [0, 1, 2].map(|i| DPoly(vec![o[i].clone(), d[i].clone()])),
394                    w: DPoly::constant(dy(1.0)),
395                    w_degree: 0,
396                })
397            }
398            Curve3::Circle(c) => Self::conic(&c.frame, c.radius, c.radius),
399            Curve3::Ellipse(e) => Self::conic(&e.frame, e.semi_axis_x, e.semi_axis_y),
400            _ => Err(ExactCurveRefusal::UnsupportedCurve),
401        }
402    }
403
404    /// `o + x a cos + y b sin` with the half-angle substitution, times
405    /// `1 + w^2`: `n = o (1 + w^2) + x a (1 - w^2) + y b (2 w)`.
406    fn conic(frame: &Frame3, a: f64, b: f64) -> Result<Self, ExactCurveRefusal> {
407        frame_finite(frame)?;
408        positive(&[a, b])?;
409        let (o, x, y) = (v3(frame.origin), v3(frame.x), v3(frame.y));
410        if cross(&x, &y).iter().all(is_zero) {
411            return Err(ExactCurveRefusal::Degenerate);
412        }
413        let (a, b) = (dy(a), dy(b));
414        let n = [0, 1, 2].map(|i| {
415            let xa = x[i].mul(&a);
416            DPoly(vec![
417                o[i].add(&xa),
418                dy(2.0).mul(&y[i]).mul(&b),
419                o[i].sub(&xa),
420            ])
421        });
422        Ok(Self {
423            n,
424            w: DPoly(vec![dy(1.0), Dyadic::zero(), dy(1.0)]),
425            w_degree: 2,
426        })
427    }
428
429    /// The homogenised linear form `u . (X - o)`, times `w`.
430    fn linear(&self, u: &V3, o: &V3) -> DPoly {
431        let mut out = self.w.scale(&dot(u, o).neg());
432        for (n, c) in self.n.iter().zip(u) {
433            out = out.add(&n.scale(c));
434        }
435        out
436    }
437}
438
439// --- the other operand as equations ----------------------------------------
440
441/// Polynomial equations `p = 0` with their degree `k` in `X`, plus an
442/// optional side condition `c >= 0` (the cone's nappe).
443struct Locus {
444    equations: Vec<(DPoly, usize)>,
445    condition: Option<DPoly>,
446}
447
448/// Local coordinates by Cramer's rule, times the frame determinant `det`:
449/// `alpha = det(P, y, z)`, `beta = det(x, P, z)`, `gamma = det(x, y, P)`.
450struct Local {
451    alpha: DPoly,
452    beta: DPoly,
453    gamma: DPoly,
454    det: Dyadic,
455}
456
457fn local(frame: (&V3, &V3, &V3, &V3), param: &Param) -> Result<Local, ExactCurveRefusal> {
458    let (o, x, y, z) = frame;
459    let det = dot(x, &cross(y, z));
460    if is_zero(&det) {
461        return Err(ExactCurveRefusal::Degenerate);
462    }
463    Ok(Local {
464        alpha: param.linear(&cross(y, z), o),
465        beta: param.linear(&cross(z, x), o),
466        gamma: param.linear(&cross(x, y), o),
467        det,
468    })
469}
470
471fn frame_of(f: &Frame3) -> Result<(V3, V3, V3, V3), ExactCurveRefusal> {
472    frame_finite(f)?;
473    Ok((v3(f.origin), v3(f.x), v3(f.y), v3(f.z)))
474}
475
476fn surface_locus(surface: &Surface, param: &Param) -> Result<Locus, ExactCurveRefusal> {
477    let one = |p: DPoly, k: usize| Locus {
478        equations: vec![(p, k)],
479        condition: None,
480    };
481    // (det * w)^2, the scale every radius term carries.
482    let dw2 = |det: &Dyadic| param.w.scale(det).square();
483    match surface {
484        Surface::Plane(p) => {
485            let (o, x, y, _) = frame_of(&p.frame)?;
486            let normal = cross(&x, &y);
487            if normal.iter().all(is_zero) {
488                return Err(ExactCurveRefusal::Degenerate);
489            }
490            Ok(one(param.linear(&normal, &o), 1))
491        }
492        Surface::Cylinder(c) => {
493            positive(&[c.radius])?;
494            let f = frame_of(&c.frame)?;
495            let l = local((&f.0, &f.1, &f.2, &f.3), param)?;
496            let r2 = dy(c.radius).square();
497            let p = l
498                .alpha
499                .square()
500                .add(&l.beta.square())
501                .sub(&dw2(&l.det).scale(&r2));
502            Ok(one(p, 2))
503        }
504        Surface::EllipticalCylinder(c) => {
505            positive(&[c.semi_axis_x, c.semi_axis_y])?;
506            let f = frame_of(&c.frame)?;
507            let l = local((&f.0, &f.1, &f.2, &f.3), param)?;
508            let (a2, b2) = (dy(c.semi_axis_x).square(), dy(c.semi_axis_y).square());
509            let p = l
510                .alpha
511                .square()
512                .scale(&b2)
513                .add(&l.beta.square().scale(&a2))
514                .sub(&dw2(&l.det).scale(&a2.mul(&b2)));
515            Ok(one(p, 2))
516        }
517        Surface::Sphere(s) => {
518            positive(&[s.radius])?;
519            let f = frame_of(&s.frame)?;
520            let l = local((&f.0, &f.1, &f.2, &f.3), param)?;
521            let r2 = dy(s.radius).square();
522            let p = l
523                .alpha
524                .square()
525                .add(&l.beta.square())
526                .add(&l.gamma.square())
527                .sub(&dw2(&l.det).scale(&r2));
528            Ok(one(p, 2))
529        }
530        Surface::Cone(c) => {
531            let slope = c.semi_angle.tan();
532            finite(&[c.radius, c.semi_angle, slope])?;
533            let f = frame_of(&c.frame)?;
534            let l = local((&f.0, &f.1, &f.2, &f.3), param)?;
535            // Local radius times det: R det w + slope gamma.
536            let reach = param
537                .w
538                .scale(&dy(c.radius).mul(&l.det))
539                .add(&l.gamma.scale(&dy(slope)));
540            let p = l.alpha.square().add(&l.beta.square()).sub(&reach.square());
541            // R + v slope >= 0 with v = gamma / (det w): the sign of
542            // `reach` times the sign of det (w is positive).
543            let condition = if l.det.sign() == Some(Sign::Negative) {
544                reach.scale(&dy(-1.0))
545            } else {
546                reach
547            };
548            Ok(Locus {
549                equations: vec![(p, 2)],
550                condition: Some(condition),
551            })
552        }
553        Surface::Torus(t) => {
554            positive(&[t.major_radius, t.minor_radius])?;
555            let f = frame_of(&t.frame)?;
556            let l = local((&f.0, &f.1, &f.2, &f.3), param)?;
557            let (big2, small2) = (dy(t.major_radius).square(), dy(t.minor_radius).square());
558            let planar = l.alpha.square().add(&l.beta.square());
559            let dw2 = dw2(&l.det);
560            // (|L|^2 + R^2 - r^2)^2 = 4 R^2 (Lx^2 + Ly^2), times det^4 w^4.
561            let s = planar
562                .add(&l.gamma.square())
563                .add(&dw2.scale(&big2.sub(&small2)));
564            let p = s.square().sub(&dw2.mul(&planar).scale(&dy(4.0).mul(&big2)));
565            Ok(one(p, 4))
566        }
567        _ => Err(ExactCurveRefusal::UnsupportedSurface),
568    }
569}
570
571fn curve_locus(curve: &Curve3, param: &Param) -> Result<Locus, ExactCurveRefusal> {
572    match curve {
573        Curve3::Line(l) => {
574            finite(&[
575                l.origin.x,
576                l.origin.y,
577                l.origin.z,
578                l.direction.x,
579                l.direction.y,
580                l.direction.z,
581            ])?;
582            let (o, d) = (v3(l.origin), v3(l.direction));
583            if d.iter().all(is_zero) {
584                return Err(ExactCurveRefusal::Degenerate);
585            }
586            // (X - o) x d = 0, one linear equation per component.
587            let z = Dyadic::zero;
588            let rows = [
589                [z(), d[2].clone(), d[1].neg()],
590                [d[2].neg(), z(), d[0].clone()],
591                [d[1].clone(), d[0].neg(), z()],
592            ];
593            Ok(Locus {
594                equations: rows.iter().map(|u| (param.linear(u, &o), 1)).collect(),
595                condition: None,
596            })
597        }
598        Curve3::Circle(c) => conic_locus(&c.frame, c.radius, c.radius, param),
599        Curve3::Ellipse(e) => conic_locus(&e.frame, e.semi_axis_x, e.semi_axis_y, param),
600        _ => Err(ExactCurveRefusal::UnsupportedCurve),
601    }
602}
603
604/// A conic is its plane (`gamma = 0` with `z = x cross y`) and, within it,
605/// `b^2 alpha^2 + a^2 beta^2 = a^2 b^2 det^2`.
606fn conic_locus(frame: &Frame3, a: f64, b: f64, param: &Param) -> Result<Locus, ExactCurveRefusal> {
607    positive(&[a, b])?;
608    let (o, x, y, _) = frame_of(frame)?;
609    let z = cross(&x, &y);
610    let l = local((&o, &x, &y, &z), param)?;
611    let (a2, b2) = (dy(a).square(), dy(b).square());
612    let ring = l
613        .alpha
614        .square()
615        .scale(&b2)
616        .add(&l.beta.square().scale(&a2))
617        .sub(&param.w.scale(&l.det).square().scale(&a2.mul(&b2)));
618    Ok(Locus {
619        equations: vec![(l.gamma, 1), (ring, 2)],
620        condition: None,
621    })
622}
623
624fn lift(curve: &Curve2) -> Result<Curve3, ExactCurveRefusal> {
625    use axiolid_curve::{Circle3, Ellipse3, Line3};
626    let frame = |f: &axiolid_core::Frame2| Frame3 {
627        origin: Point3::new(f.origin.x, f.origin.y, 0.0),
628        x: Vec3::new(f.x.x, f.x.y, 0.0),
629        y: Vec3::new(f.y.x, f.y.y, 0.0),
630        z: Vec3::Z,
631    };
632    Ok(match curve {
633        Curve2::Line(l) => Curve3::Line(Line3 {
634            origin: Point3::new(l.origin.x, l.origin.y, 0.0),
635            direction: Vec3::new(l.direction.x, l.direction.y, 0.0),
636        }),
637        Curve2::Circle(c) => Curve3::Circle(Circle3 {
638            frame: frame(&c.frame),
639            radius: c.radius,
640        }),
641        Curve2::Ellipse(e) => Curve3::Ellipse(Ellipse3 {
642            frame: frame(&e.frame),
643            semi_axis_x: e.semi_axis_x,
644            semi_axis_y: e.semi_axis_y,
645        }),
646        _ => return Err(ExactCurveRefusal::UnsupportedCurve),
647    })
648}
649
650// --- solving -----------------------------------------------------------------
651
652/// Multiplicity of `root` in `g`: the order of the first derivative that
653/// does not vanish there.
654fn multiplicity(root: &RealRoot, g: &IntPoly) -> usize {
655    let mut order = 1;
656    let mut q = g.derivative();
657    while !q.is_zero() && root.sign_of(&q) == Sign::Zero {
658        order += 1;
659        q = q.derivative();
660    }
661    order
662}
663
664fn solve(curve: &Curve3, param: &Param, locus: &Locus) -> ExactCurveIntersection {
665    // Each equation becomes an integer polynomial whose FORMAL degree is
666    // k * deg(w): its coefficient of s^(k deg w) is the equation at the
667    // conic's point at infinity (theta = pi).
668    let polys: Vec<(IntPoly, usize)> = locus
669        .equations
670        .iter()
671        .map(|(p, k)| (p.to_int(), k * param.w_degree))
672        .filter(|(p, _)| !p.is_zero())
673        .collect();
674    if polys.is_empty() {
675        return ExactCurveIntersection::Contained;
676    }
677    // Common roots of all equations are the roots of their gcd.
678    let g = polys
679        .iter()
680        .skip(1)
681        .fold(polys[0].0.clone(), |acc, (p, _)| acc.gcd(p));
682    let condition_int = locus.condition.as_ref().map(DPoly::to_int);
683    let allowed = |root: &RealRoot| {
684        condition_int
685            .as_ref()
686            .is_none_or(|c| c.is_zero() || root.sign_of(c) != Sign::Negative)
687    };
688
689    let is_line = param.w_degree == 0;
690    let mut hits: Vec<(i8, ExactCurveHit)> = Vec::new();
691    for root in g.real_roots() {
692        if !allowed(&root) {
693            continue;
694        }
695        let multiplicity = multiplicity(&root, &g);
696        // Order along the conic by theta in [0, 2 pi): w >= 0 first, then
697        // the antipode, then w < 0.
698        let band = if is_line || root.cmp_dyadic(&Dyadic::zero()) != Sign::Negative {
699            0
700        } else {
701            2
702        };
703        let parameter = if is_line {
704            ExactCurveParameter::Line(root)
705        } else {
706            ExactCurveParameter::HalfAngle(root)
707        };
708        hits.push((band, hit(curve, parameter, multiplicity)));
709    }
710    if !is_line {
711        // A common zero at infinity: every equation lost degree. Its
712        // multiplicity is the smallest loss.
713        let lost = polys
714            .iter()
715            .map(|(p, formal)| formal - p.degree().unwrap_or(0))
716            .min()
717            .unwrap_or(0);
718        let on_nappe = locus.condition.as_ref().is_none_or(|c| {
719            // The condition's formal degree is deg(w) = 2; its leading
720            // coefficient is its value at the antipode.
721            c.0.get(2)
722                .is_none_or(|lead| lead.sign() != Some(Sign::Negative))
723        });
724        if lost >= 1 && on_nappe {
725            hits.push((1, hit(curve, ExactCurveParameter::Antipode, lost)));
726        }
727    }
728    // Roots come in increasing order; a stable sort by band keeps it.
729    hits.sort_by_key(|(band, _)| *band);
730    ExactCurveIntersection::Points(hits.into_iter().map(|(_, h)| h).collect())
731}
732
733fn hit(curve: &Curve3, parameter: ExactCurveParameter, multiplicity: usize) -> ExactCurveHit {
734    let point =
735        evaluate3(curve, parameter.approx()).unwrap_or(Point3::new(f64::NAN, f64::NAN, f64::NAN));
736    ExactCurveHit {
737        parameter,
738        multiplicity,
739        point,
740    }
741}
axiolid_nurbs/exact_curve_intersection.rs

axiolid_nurbs/
exact_curve_intersection.rs
