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

axiolid_nurbs/
exact_curve_intersection.rs
