axiolid_nurbs/
implicit_ops.rs

1//! Working with section curves through their implicit form (ADR 0077):
2//! pcurves of a section on any analytic face, and a curve of the section
3//! families against a surface.
4//!
5//! Every section family is a field's zero set on a carrier: an
6//! `ImplicitSection3` is one by definition. A ruled section (ADR 0076) is
7//! the root branch of `a v^2 + b v + c = 0` on its ruled carrier, and a
8//! torus section is where `A cos u + B sin u = C` on its torus. Each is
9//! traced once in its own carrier's parameters. The stretch that matches
10//! the curve is cut out by its end points, so the operations below have one
11//! implementation for all of them.
12
13use axiolid_core::{Interval, Point2, Scalar, Vec2};
14use axiolid_curve::implicit::{bound_simple, partial, Cell, Range};
15use axiolid_curve::{
16    Axis, Basis, Carrier, Curve3, Field2, ImplicitCell, ImplicitCurve2, ImplicitSection3,
17    SeriesField2, Trig2,
18};
19use axiolid_surface::Surface;
20use core::f64::consts::{PI, TAU};
21
22use crate::exact_curve_intersection::{
23    ExactCurveHit, ExactCurveIntersection, ExactCurveParameter, ExactCurveRefusal, Isolated,
24};
25use crate::exact_surface_intersection::ExactIntersectionRefusal;
26use crate::field::{carrier_of, section_field};
27use crate::implicit_trace::Periodic;
28
29/// Every component of the section of `surface` by `other` in a window of
30/// `surface`'s parameters, as curves in those parameters.
31///
32/// A window a whole turn wide in a periodic parameter wraps.
33///
34/// # Errors
35///
36/// As [`crate::implicit_surface_intersection`]; `Disjoint` is not an
37/// error here (the list is empty).
38pub fn trace_section_pcurves(
39    surface: &Surface,
40    other: &Surface,
41    window: (Point2, Point2),
42) -> Result<Vec<ImplicitCurve2>, ExactIntersectionRefusal> {
43    let carrier = carrier_of(surface).ok_or(ExactIntersectionRefusal::UnsupportedPair)?;
44    let field = section_field(&carrier, other).ok_or(ExactIntersectionRefusal::UnsupportedPair)?;
45    // A quarter of a section's budget: the section itself is traced
46    // already, and a face without these pcurves reads it back instead.
47    trace_field(&carrier, &field, window, crate::implicit_trace::BUDGET / 4)
48}
49
50fn trace_field(
51    carrier: &Carrier,
52    field: &Field2,
53    window: (Point2, Point2),
54    budget: usize,
55) -> Result<Vec<ImplicitCurve2>, ExactIntersectionRefusal> {
56    let (pu, pv) = carrier.periodic();
57    let (lo, hi) = window;
58    let wraps = |periodic: bool, a: Scalar, b: Scalar| {
59        periodic && ((b - a) - TAU).abs() <= 1e-12 * (1.0 + a.abs())
60    };
61    let periodic = Periodic {
62        u: wraps(pu, lo.x, hi.x),
63        v: wraps(pv, lo.y, hi.y),
64    };
65    crate::implicit_trace::trace_within(field, Cell { lo, hi }, periodic, budget)
66        .map_err(crate::implicit_section::refusal_of)
67}
68
69/// The stretch of one of `curves` that runs from `start` through `first`
70/// and then `second` to `end` (all in parameters, each possibly a whole
71/// period off), as a curve of its own running from `start` to `end`.
72/// `closed` asks for a whole loop from `start` back to itself, in the
73/// direction that meets `first` before `second`.
74#[must_use]
75pub fn extract_stretch(
76    curves: &[ImplicitCurve2],
77    periodic: (bool, bool),
78    start: Point2,
79    [first, second]: [Point2; 2],
80    end: Point2,
81    closed: bool,
82) -> Option<ImplicitCurve2> {
83    let shifts = |p: Point2| -> Vec<Point2> {
84        let ku: &[Scalar] = if periodic.0 {
85            &[0.0, -1.0, 1.0, -2.0, 2.0]
86        } else {
87            &[0.0]
88        };
89        let kv: &[Scalar] = if periodic.1 {
90            &[0.0, -1.0, 1.0, -2.0, 2.0]
91        } else {
92            &[0.0]
93        };
94        let mut out = Vec::new();
95        for a in ku {
96            for b in kv {
97                out.push(Point2::new(p.x + a * TAU, p.y + b * TAU));
98            }
99        }
100        out
101    };
102    let find = |curve: &ImplicitCurve2, p: Point2| -> Option<Scalar> {
103        let eps = 1e-7 * (1.0 + p.x.abs().max(p.y.abs()));
104        let mut best: Option<(Scalar, Scalar)> = None;
105        for q in shifts(p) {
106            let Some(t) = curve.parameter_of(q) else {
107                continue;
108            };
109            let Some(on) = curve.point(t) else {
110                continue;
111            };
112            // Compare up to whole periods.
113            let mut d = on - q;
114            if periodic.0 {
115                d.x -= (d.x / TAU).round() * TAU;
116            }
117            if periodic.1 {
118                d.y -= (d.y / TAU).round() * TAU;
119            }
120            let miss = d.length();
121            if miss <= eps && best.is_none_or(|(m, _)| miss < m) {
122                best = Some((miss, t));
123            }
124        }
125        best.map(|(_, t)| t)
126    };
127    for curve in curves {
128        let (Some(ts), Some(t1), Some(t2), Some(te)) = (
129            find(curve, start),
130            find(curve, first),
131            find(curve, second),
132            find(curve, end),
133        ) else {
134            continue;
135        };
136        let n = curve.end();
137        let closure = curve.closure(periodic.0, periodic.1);
138        let ahead = |t: Scalar| (t - ts).rem_euclid(n);
139        let stretch = if closed {
140            // The whole loop from `start` round to itself.
141            let whole = curve.rotated(ts, closure?)?;
142            if ahead(t1) < ahead(t2) {
143                whole
144            } else {
145                whole.reversed()
146            }
147        } else if let Some(closure) = closure {
148            // On a loop: forward if the inner points come in order before
149            // `end` going forward from `start`.
150            if ahead(t1) < ahead(t2) && ahead(t2) < ahead(te) {
151                curve.sub(ts, te, Some(closure))?
152            } else {
153                curve.sub(te, ts, Some(closure))?.reversed()
154            }
155        } else if ts <= t1 && t1 <= t2 && t2 <= te {
156            curve.sub(ts, te, None)?
157        } else if te <= t2 && t2 <= t1 && t1 <= ts {
158            curve.sub(te, ts, None)?.reversed()
159        } else {
160            continue;
161        };
162        return Some(stretch);
163    }
164    None
165}
166
167/// The carrier and defining field of a section family, or `None` for
168/// other curves.
169fn defining_field(curve: &Curve3) -> Option<(Carrier, Field2)> {
170    let fourier = |t: &Trig2| vec![t.constant, t.cos, t.sin, t.cos2, t.sin2];
171    match curve {
172        Curve3::ImplicitSection(s) => Some((s.carrier.clone(), s.curve.field.clone())),
173        Curve3::RuledSection(r) => {
174            // c(u) + b(u) v + a(u) v^2, u harmonic, v a power.
175            let (a, b, c) = (
176                fourier(&r.graph.a),
177                fourier(&r.graph.b),
178                fourier(&r.graph.c),
179            );
180            let coefficients = (0..5).map(|i| vec![c[i], b[i], a[i]]).collect();
181            Some((
182                Carrier::Ruled(r.carrier),
183                Field2::Series(SeriesField2 {
184                    u: Basis::Fourier,
185                    v: Basis::Power,
186                    coefficients,
187                }),
188            ))
189        }
190        Curve3::TorusSection(t) => {
191            // A(v) cos u + B(v) sin u - C(v), both harmonic.
192            let (a, b, c) = (
193                fourier(&t.graph.a),
194                fourier(&t.graph.b),
195                fourier(&t.graph.c),
196            );
197            let coefficients = vec![c.iter().map(|x| -x).collect(), a, b];
198            Some((
199                Carrier::Torus(t.torus),
200                Field2::Series(SeriesField2 {
201                    u: Basis::Fourier,
202                    v: Basis::Fourier,
203                    coefficients,
204                }),
205            ))
206        }
207        _ => None,
208    }
209}
210
211/// A section-family curve over `span` as an [`ImplicitSection3`], with the
212/// same point set and direction; `None` for other families or where the
213/// trace is refused.
214#[must_use]
215pub fn implicit_view(curve: &Curve3, span: Interval) -> Option<ImplicitSection3> {
216    if let Curve3::ImplicitSection(s) = curve {
217        let closed = (span.end - span.start).abs() >= s.curve.end() - 1e-12;
218        let sub = if closed {
219            s.curve.clone()
220        } else {
221            let (a, b) = (span.start.min(span.end), span.start.max(span.end));
222            let piece = s.curve.sub(a, b, None)?;
223            if span.end < span.start {
224                piece.reversed()
225            } else {
226                piece
227            }
228        };
229        return Some(ImplicitSection3 {
230            carrier: s.carrier.clone(),
231            curve: sub,
232        });
233    }
234    let (carrier, field) = defining_field(curve)?;
235    let at = |t: Scalar| -> Option<Point2> {
236        let p = axiolid_evaluate::evaluate3(curve, t).ok()?;
237        let (u, v) = carrier.parameters(p);
238        Some(Point2::new(u, v))
239    };
240    // Carrier parameters along the span, unwrapped, to size the window.
241    let n = 64;
242    let mut samples = Vec::with_capacity(n + 1);
243    for i in 0..=n {
244        let t = span.start + (span.end - span.start) * i as Scalar / n as Scalar;
245        let mut p = at(t)?;
246        if let Some(q) = samples.last().copied() {
247            let q: Point2 = q;
248            let (pu, pv) = carrier.periodic();
249            if pu {
250                p.x += ((q.x - p.x) / TAU).round() * TAU;
251            }
252            if pv {
253                p.y += ((q.y - p.y) / TAU).round() * TAU;
254            }
255        }
256        samples.push(p);
257    }
258    let (mut lo, mut hi) = (samples[0], samples[0]);
259    for p in &samples {
260        lo = lo.min(*p);
261        hi = hi.max(*p);
262    }
263    let pad = (hi - lo) * 0.25 + Vec2::splat(0.05);
264    let (mut lo, mut hi) = (lo - pad, hi + pad);
265    // Never more than a turn in a periodic parameter.
266    let (pu, pv) = carrier.periodic();
267    if pu && hi.x - lo.x >= TAU {
268        let c = 0.5 * (lo.x + hi.x);
269        (lo.x, hi.x) = (c - PI, c + PI);
270    }
271    if pv && hi.y - lo.y >= TAU {
272        let c = 0.5 * (lo.y + hi.y);
273        (lo.y, hi.y) = (c - PI, c + PI);
274    }
275    let curves = trace_field(&carrier, &field, (lo, hi), crate::implicit_trace::BUDGET).ok()?;
276    // A loop closes in space, whatever whole turns its parameters make.
277    let (first, last) = (
278        axiolid_evaluate::evaluate3(curve, span.start).ok()?,
279        axiolid_evaluate::evaluate3(curve, span.end).ok()?,
280    );
281    let closed = (first - last).length() <= 1e-9 * (1.0 + first.length());
282    let stretch = extract_stretch(
283        &curves,
284        (pu, pv),
285        samples[0],
286        [samples[n / 3], samples[2 * n / 3]],
287        samples[n],
288        closed,
289    )?;
290    Some(ImplicitSection3 {
291        carrier,
292        curve: stretch,
293    })
294}
295
296/// Where a section-family curve meets a surface, with parameters on the
297/// curve itself.
298///
299/// The curve is read as an implicit curve on its carrier (over `span`, or
300/// its whole traced extent for an [`ImplicitSection3`]); the surface's
301/// equation on that carrier is a second field, and its roots along each
302/// cell are isolated with interval bounds over boxes certain to hold the
303/// curve. A root where the second field changes sign is a crossing; one
304/// where it only touches zero is reported with multiplicity 2; a curve on
305/// which the second field vanishes identically is `Contained`.
306///
307/// # Errors
308///
309/// `UnsupportedCurve` for a curve that is not a section family or whose
310/// trace was refused, `UnsupportedSurface` for a B-spline surface.
311pub fn section_curve_surface_intersection(
312    curve: &Curve3,
313    span: Interval,
314    surface: &Surface,
315) -> Result<ExactCurveIntersection, ExactCurveRefusal> {
316    let view = implicit_view(curve, span).ok_or(ExactCurveRefusal::UnsupportedCurve)?;
317    let other =
318        section_field(&view.carrier, surface).ok_or(ExactCurveRefusal::UnsupportedSurface)?;
319    let own = &view.curve;
320    // Contained: the second field vanishes along the whole curve, to the
321    // rounding of its own size.
322    let scale = 1e-9 * other.magnitude().max(1.0);
323    let samples: Vec<Scalar> = (0..=32)
324        .filter_map(|i| {
325            let p = own.point(own.end() * i as Scalar / 32.0)?;
326            Some(other.value(p))
327        })
328        .collect();
329    if samples.len() == 33 && samples.iter().all(|h| h.abs() <= scale) {
330        return Ok(ExactCurveIntersection::Contained);
331    }
332    let mut roots = roots_along(own, &other);
333    roots.sort_by(|a, b| a.0.total_cmp(&b.0));
334    roots.dedup_by(|a, b| (a.0 - b.0).abs() <= 1e-9);
335    let mut hits = Vec::new();
336    for (t, multiplicity) in roots {
337        let Some(p) = own.point(t) else { continue };
338        let point = view.carrier.jet(p.x, p.y).point;
339        // Back to the curve's own parameter.
340        let parameter = match curve {
341            Curve3::ImplicitSection(_) => {
342                if span.end < span.start {
343                    span.start - t
344                } else {
345                    span.start + t
346                }
347            }
348            _ => {
349                match axiolid_evaluate::curve::invert3(curve, point, axiolid_core::Tolerance::METRE)
350                {
351                    Ok(x) => nearest_turn(x, span),
352                    Err(_) => continue,
353                }
354            }
355        };
356        hits.push(ExactCurveHit {
357            parameter: ExactCurveParameter::Certified(Isolated::new(parameter)),
358            multiplicity,
359            point,
360        });
361    }
362    Ok(ExactCurveIntersection::Points(hits))
363}
364
365/// The roots of `field` along an implicit curve, by the curve's parameter,
366/// with multiplicity 1 where the field changes sign and 2 where it only
367/// touches zero. Each cell's stretch is bounded by a box certain to hold it
368/// (the solved parameter moves at most `max |F_free| / min |F_solved|` per
369/// unit of the free one), and the field's roots are isolated with interval
370/// bounds over those boxes.
371pub(crate) fn roots_along(own: &ImplicitCurve2, other: &Field2) -> Vec<(Scalar, usize)> {
372    let d_u = partial(&own.field, true);
373    let d_v = partial(&own.field, false);
374    let h_u = partial(other, true);
375    let h_v = partial(other, false);
376    let mut roots: Vec<(Scalar, usize)> = Vec::new();
377    for (index, cell) in own.cells.iter().enumerate() {
378        let job = CellRoots {
379            curve: own,
380            cell,
381            other,
382            d_u: &d_u,
383            d_v: &d_v,
384            h_u: &h_u,
385            h_v: &h_v,
386        };
387        let mut found = Vec::new();
388        job.roots(0.0, 1.0, 0, &mut found);
389        for (s, m) in found {
390            roots.push((index as Scalar + s, m));
391        }
392    }
393    roots.sort_by(|a, b| a.0.total_cmp(&b.0));
394    roots.dedup_by(|a, b| (a.0 - b.0).abs() <= 1e-9);
395    roots
396}
397
398/// `x` moved by whole turns into `span` when that is possible.
399fn nearest_turn(x: Scalar, span: Interval) -> Scalar {
400    let (lo, hi) = (span.start.min(span.end), span.start.max(span.end));
401    for k in [0.0, 1.0, -1.0, 2.0, -2.0] {
402        let y = x + k * TAU;
403        if y >= lo - 1e-9 && y <= hi + 1e-9 {
404            return y;
405        }
406    }
407    x
408}
409
410/// The roots of a second field along one cell of an implicit curve.
411struct CellRoots<'a> {
412    curve: &'a ImplicitCurve2,
413    cell: &'a ImplicitCell,
414    other: &'a Field2,
415    d_u: &'a Field2,
416    d_v: &'a Field2,
417    h_u: &'a Field2,
418    h_v: &'a Field2,
419}
420
421impl CellRoots<'_> {
422    fn free(&self, s: Scalar) -> Scalar {
423        self.cell.from + (self.cell.to - self.cell.from) * s
424    }
425
426    fn place(&self, free: Scalar, solved: Scalar) -> Point2 {
427        match self.cell.axis {
428            Axis::U => Point2::new(free, solved),
429            Axis::V => Point2::new(solved, free),
430        }
431    }
432
433    fn solved(&self, free: Scalar) -> Option<Scalar> {
434        if self.cell.bridge.is_some() {
435            return self.curve.solve_cell(self.cell, free);
436        }
437        let one = ImplicitCurve2 {
438            field: self.curve.field.clone(),
439            cells: vec![ImplicitCell {
440                from: free,
441                to: free,
442                ..*self.cell
443            }],
444        };
445        let p = one.point(0.0)?;
446        Some(match self.cell.axis {
447            Axis::U => p.y,
448            Axis::V => p.x,
449        })
450    }
451
452    fn h(&self, s: Scalar) -> Option<Scalar> {
453        let free = self.free(s);
454        Some(self.other.value(self.place(free, self.solved(free)?)))
455    }
456
457    /// A box certain to hold the curve for `s` in `[s0, s1]`, and bounds of
458    /// the solved parameter's slope there.
459    fn hull(&self, s0: Scalar, s1: Scalar) -> Option<(Cell, Range)> {
460        let (f0, f1) = (self.free(s0), self.free(s1));
461        let (f_lo, f_hi) = (f0.min(f1), f0.max(f1));
462        let w0 = self.solved(f0)?;
463        let make = |w_lo: Scalar, w_hi: Scalar| {
464            let (a, b) = (self.place(f_lo, w_lo), self.place(f_hi, w_hi));
465            Cell {
466                lo: a.min(b),
467                hi: a.max(b),
468            }
469        };
470        let (d_free, d_solved) = match self.cell.axis {
471            Axis::U => (self.d_u, self.d_v),
472            Axis::V => (self.d_v, self.d_u),
473        };
474        let whole = make(self.cell.low, self.cell.high);
475        let free = bound_simple(d_free, &whole);
476        let solved = bound_simple(d_solved, &whole);
477        if solved.straddles_zero() {
478            return None;
479        }
480        let floor = solved.lo.abs().min(solved.hi.abs());
481        let big = free.lo.abs().max(free.hi.abs());
482        let reach = big / floor * (f_hi - f_lo);
483        let hull = make(
484            (w0 - reach).max(self.cell.low),
485            (w0 + reach).min(self.cell.high),
486        );
487        // Slope bound of the solved parameter over the hull: -F_free / F_solved.
488        let free = bound_simple(d_free, &hull);
489        let solved = bound_simple(d_solved, &hull);
490        let q = [
491            -free.lo / solved.lo,
492            -free.lo / solved.hi,
493            -free.hi / solved.lo,
494            -free.hi / solved.hi,
495        ];
496        let slope = Range {
497            lo: q.iter().copied().fold(Scalar::INFINITY, Scalar::min),
498            hi: q.iter().copied().fold(Scalar::NEG_INFINITY, Scalar::max),
499        };
500        Some((hull, slope))
501    }
502
503    fn roots(&self, s0: Scalar, s1: Scalar, depth: u32, out: &mut Vec<(Scalar, usize)>) {
504        // A bridge into a crossing is too short and too close to the
505        // singular point for bounds: its sign changes are scanned.
506        if self.cell.bridge.is_some() {
507            let n = 64;
508            let at = |k: usize| s0 + (s1 - s0) * k as Scalar / n as Scalar;
509            let mut last = self.h(at(0));
510            for k in 1..=n {
511                let now = self.h(at(k));
512                if let (Some(a), Some(b)) = (last, now) {
513                    if a == 0.0 && k == 1 {
514                        out.push((at(0), 1));
515                    } else if (a < 0.0) != (b < 0.0) && b != 0.0 {
516                        out.push((self.bisect(at(k - 1), at(k), a), 1));
517                    } else if b == 0.0 {
518                        out.push((at(k), 1));
519                    }
520                }
521                last = now;
522            }
523            return;
524        }
525        let Some((hull, slope)) = self.hull(s0, s1) else {
526            return;
527        };
528        if !bound_simple(self.other, &hull).straddles_zero() {
529            return;
530        }
531        let (Some(a), Some(b)) = (self.h(s0), self.h(s1)) else {
532            return;
533        };
534        // d h / d free = H_free + H_solved * slope.
535        let (h_free, h_solved) = match self.cell.axis {
536            Axis::U => (self.h_u, self.h_v),
537            Axis::V => (self.h_v, self.h_u),
538        };
539        let hf = bound_simple(h_free, &hull);
540        let hs = bound_simple(h_solved, &hull);
541        let p = [
542            hs.lo * slope.lo,
543            hs.lo * slope.hi,
544            hs.hi * slope.lo,
545            hs.hi * slope.hi,
546        ];
547        let lo = hf.lo + p.iter().copied().fold(Scalar::INFINITY, Scalar::min);
548        let hi = hf.hi + p.iter().copied().fold(Scalar::NEG_INFINITY, Scalar::max);
549        let monotone = lo > 0.0 || hi < 0.0;
550        if monotone {
551            if a == 0.0 {
552                out.push((s0, 1));
553            } else if (a < 0.0) != (b < 0.0) {
554                out.push((self.bisect(s0, s1, a), 1));
555            }
556            return;
557        }
558        if depth >= 44 || s1 - s0 <= 1e-13 {
559            if (a < 0.0) != (b < 0.0) {
560                out.push((self.bisect(s0, s1, a), 1));
561            } else {
562                // Touching: the second field reaches zero without crossing.
563                let scale = 1e-9 * self.other.magnitude().max(1.0);
564                if a.abs().min(b.abs()) <= scale {
565                    out.push((0.5 * (s0 + s1), 2));
566                }
567            }
568            return;
569        }
570        let m = 0.5 * (s0 + s1);
571        self.roots(s0, m, depth + 1, out);
572        self.roots(m, s1, depth + 1, out);
573    }
574
575    fn bisect(&self, mut s0: Scalar, mut s1: Scalar, a: Scalar) -> Scalar {
576        let negative = a < 0.0;
577        for _ in 0..100 {
578            let m = 0.5 * (s0 + s1);
579            if m <= s0 || m >= s1 {
580                break;
581            }
582            match self.h(m) {
583                Some(0.0) => return m,
584                Some(h) if (h < 0.0) == negative => s0 = m,
585                Some(_) => s1 = m,
586                None => break,
587            }
588        }
589        0.5 * (s0 + s1)
590    }
591}
592
593/// A line, circle or ellipse against a B-spline surface: the curve is the
594/// meeting of two analytic surfaces (a line of two planes, a circle of its
595/// plane and a sphere, an ellipse of its plane and an elliptical cylinder).
596/// The first is traced on the spline (certified, over its Bernstein
597/// patches); the second's roots along each traced piece are the hits.
598///
599/// # Errors
600///
601/// `UnsupportedCurve` for other families, `UnsupportedSurface` for a
602/// surface that is not a B-spline or a spline whose trace is refused.
603pub fn conic_spline_intersection(
604    curve: &Curve3,
605    surface: &Surface,
606) -> Result<ExactCurveIntersection, ExactCurveRefusal> {
607    use axiolid_core::Frame3;
608    use axiolid_surface::{EllipticalCylinder, Plane, Sphere};
609    if !matches!(surface, Surface::BSpline(_)) {
610        return Err(ExactCurveRefusal::UnsupportedSurface);
611    }
612    let plane = |origin: axiolid_core::Point3, z: axiolid_core::Vec3| -> Surface {
613        let z = z.normalize();
614        let helper = if z.x.abs() < 0.9 {
615            axiolid_core::Vec3::X
616        } else {
617            axiolid_core::Vec3::Y
618        };
619        let x = helper.cross(z).normalize();
620        Surface::Plane(Plane {
621            frame: Frame3 {
622                origin,
623                x,
624                y: z.cross(x),
625                z,
626            },
627        })
628    };
629    let (first, second) = match curve {
630        Curve3::Line(l) => {
631            let d = l.direction.normalize();
632            let helper = if d.x.abs() < 0.9 {
633                axiolid_core::Vec3::X
634            } else {
635                axiolid_core::Vec3::Y
636            };
637            let n1 = d.cross(helper).normalize();
638            (plane(l.origin, n1), plane(l.origin, d.cross(n1)))
639        }
640        Curve3::Circle(c) => (
641            plane(c.frame.origin, c.frame.x.cross(c.frame.y)),
642            Surface::Sphere(Sphere {
643                frame: c.frame,
644                radius: c.radius,
645            }),
646        ),
647        Curve3::Ellipse(e) => (
648            plane(e.frame.origin, e.frame.x.cross(e.frame.y)),
649            Surface::EllipticalCylinder(EllipticalCylinder {
650                frame: Frame3 {
651                    z: e.frame.x.cross(e.frame.y).normalize(),
652                    ..e.frame
653                },
654                semi_axis_x: e.semi_axis_x,
655                semi_axis_y: e.semi_axis_y,
656            }),
657        ),
658        _ => return Err(ExactCurveRefusal::UnsupportedCurve),
659    };
660    let pieces = match crate::implicit_section::implicit_surface_intersection(surface, &first, None)
661    {
662        Ok(pieces) => pieces,
663        Err(ExactIntersectionRefusal::Disjoint) => {
664            return Ok(ExactCurveIntersection::Points(Vec::new()))
665        }
666        Err(_) => return Err(ExactCurveRefusal::UnsupportedSurface),
667    };
668    let mut hits = Vec::new();
669    for piece in pieces {
670        let end = piece.curve.end();
671        let traced = Curve3::ImplicitSection(piece);
672        match section_curve_surface_intersection(&traced, Interval::new(0.0, end), &second)? {
673            ExactCurveIntersection::Points(found) => {
674                for hit in found {
675                    let t = match curve {
676                        Curve3::Line(l) => {
677                            (hit.point - l.origin).dot(l.direction) / l.direction.length_squared()
678                        }
679                        _ => match axiolid_evaluate::curve::invert3(
680                            curve,
681                            hit.point,
682                            axiolid_core::Tolerance::METRE,
683                        ) {
684                            Ok(t) => t,
685                            Err(_) => continue,
686                        },
687                    };
688                    hits.push(ExactCurveHit {
689                        parameter: ExactCurveParameter::Certified(Isolated::new(t)),
690                        ..hit
691                    });
692                }
693            }
694            // The traced piece lies on the second surface too: the curve
695            // runs in the spline.
696            _ => return Ok(ExactCurveIntersection::Contained),
697        }
698    }
699    hits.sort_by(|a, b| a.parameter.approx().total_cmp(&b.parameter.approx()));
700    Ok(ExactCurveIntersection::Points(hits))
701}
axiolid_nurbs/implicit_ops.rs

axiolid_nurbs/
implicit_ops.rs
