axiolid_nurbs/
section_curves.rs

1//! Curve/curve intersection for the section families (#119, B8;
2//! ADR 0077).
3//!
4//! **In the plane.** Every plane curve of the section families is a
5//! field's zero set: a line or conic by its implicit equation, a
6//! [`Sinusoid2`] by `v - w(u)`, a [`QuadraticGraph2`] by its quadratic, an
7//! [`AngleGraph2`] by `A cos u + B sin u - C`, an [`ImplicitCurve2`] by its
8//! own field. One piece is taken as traced cells, over its span (traced
9//! with certified topology when it is not already an implicit curve), and
10//! the other's field has its roots isolated along those cells with interval
11//! bounds. Each root is kept when it lies on both pieces within their
12//! spans: a graph's other branch shares its field, and is filtered here.
13//!
14//! **In space.** Two space curves meet only where one crosses a surface the
15//! other lies on: the first is intersected with the second's carrier (a
16//! traced section's surface, a conic's plane, a line's two planes) by the
17//! certified curve/surface routines, and each point is kept when it lies on
18//! the second curve within `tolerance`. In space the question is itself
19//! only meaningful up to a tolerance: two curves built from rounded
20//! numbers rarely meet exactly.
21//!
22//! [`Sinusoid2`]: axiolid_curve::Sinusoid2
23//! [`QuadraticGraph2`]: axiolid_curve::QuadraticGraph2
24//! [`AngleGraph2`]: axiolid_curve::AngleGraph2
25//! [`ImplicitCurve2`]: axiolid_curve::ImplicitCurve2
26
27use axiolid_core::{Frame3, Interval, Point2, Point3, Scalar, Tolerance, Vec2, Vec3};
28use axiolid_curve::{Basis, Carrier, Curve2, Curve3, Field2, ImplicitCurve2, SeriesField2, Trig2};
29use axiolid_evaluate::curve::{evaluate2, evaluate3, locate2, locate3};
30use axiolid_surface::{Plane, Surface};
31
32use crate::exact_curve_intersection::{ExactCurveIntersection, ExactCurveRefusal};
33
34/// One point where two curves meet, with its parameter on each.
35#[derive(Debug, Clone, Copy, PartialEq)]
36pub struct CurveCurveHit {
37    /// Parameter on the first curve.
38    pub first: Scalar,
39    /// Parameter on the second curve.
40    pub second: Scalar,
41    /// The point (on the first curve at `first`).
42    pub point: Point3,
43    /// 1 where the curves cross, 2 or more where they touch (the order of
44    /// the contact, up to eight, for a traced curve). At an end of the
45    /// traced curve, where the curve may stop short of a contact, it is the
46    /// order the enclosures see there, not a class: a touch a rounding's
47    /// width from the end can read as a crossing.
48    pub multiplicity: usize,
49}
50
51fn series(u: Basis, v: Basis, coefficients: Vec<Vec<Scalar>>) -> Field2 {
52    Field2::Series(SeriesField2 { u, v, coefficients })
53}
54
55fn fourier(t: &Trig2) -> Vec<Scalar> {
56    vec![t.constant, t.cos, t.sin, t.cos2, t.sin2]
57}
58
59/// A plane curve's defining field, or `None` for a family without one.
60fn plane_field(curve: &Curve2) -> Option<Field2> {
61    let p = Basis::Power;
62    let f = Basis::Fourier;
63    Some(match curve {
64        Curve2::Line(l) => {
65            let n = Vec2::new(-l.direction.y, l.direction.x);
66            if n.length() == 0.0 {
67                return None;
68            }
69            let n = n.normalize();
70            series(p, p, vec![vec![-n.dot(l.origin), n.y], vec![n.x, 0.0]])
71        }
72        Curve2::Circle(c) => conic_field(c.frame, c.radius, c.radius)?,
73        Curve2::Ellipse(e) => conic_field(e.frame, e.semi_axis_x, e.semi_axis_y)?,
74        Curve2::Sinusoid(w) => series(
75            f,
76            p,
77            vec![vec![-w.mean, 1.0], vec![-w.cosine, 0.0], vec![-w.sine, 0.0]],
78        ),
79        Curve2::QuadraticGraph(g) => {
80            let (a, b, c) = (fourier(&g.a), fourier(&g.b), fourier(&g.c));
81            series(f, p, (0..5).map(|i| vec![c[i], b[i], a[i]]).collect())
82        }
83        Curve2::AngleGraph(g) => {
84            let (a, b, c) = (fourier(&g.a), fourier(&g.b), fourier(&g.c));
85            series(f, f, vec![c.iter().map(|x| -x).collect(), a, b])
86        }
87        Curve2::Implicit(c) => c.field.clone(),
88        _ => return None,
89    })
90}
91
92/// `b^2 x^2 + a^2 y^2 - a^2 b^2` in an orthonormal frame's coordinates.
93fn conic_field(frame: axiolid_core::Frame2, a: Scalar, b: Scalar) -> Option<Field2> {
94    let (x, y) = (frame.x, frame.y);
95    if (x.length() - 1.0).abs() > 1e-12
96        || (y.length() - 1.0).abs() > 1e-12
97        || x.dot(y).abs() > 1e-12
98    {
99        return None;
100    }
101    // x_l = X . (p - o), y_l = Y . (p - o): affine in (u, v).
102    let lin = |axis: Vec2| [-axis.dot(frame.origin), axis.x, axis.y];
103    let (lx, ly) = (lin(x), lin(y));
104    // Square of c0 + c1 u + c2 v, as coefficients [i][j] of u^i v^j.
105    let square = |l: [Scalar; 3]| {
106        let mut c = vec![vec![0.0; 3]; 3];
107        c[0][0] = l[0] * l[0];
108        c[1][0] = 2.0 * l[0] * l[1];
109        c[0][1] = 2.0 * l[0] * l[2];
110        c[2][0] = l[1] * l[1];
111        c[1][1] = 2.0 * l[1] * l[2];
112        c[0][2] = l[2] * l[2];
113        c
114    };
115    let (sx, sy) = (square(lx), square(ly));
116    let mut c = vec![vec![0.0; 3]; 3];
117    for i in 0..3 {
118        for j in 0..3 {
119            c[i][j] = b * b * sx[i][j] + a * a * sy[i][j];
120        }
121    }
122    c[0][0] -= a * a * b * b;
123    Some(series(Basis::Power, Basis::Power, c))
124}
125
126/// A plane curve's piece as traced cells of its own field, over a window
127/// holding it.
128fn traced_piece(curve: &Curve2, span: Interval) -> Option<ImplicitCurve2> {
129    if let Curve2::Implicit(c) = curve {
130        let (a, b) = (span.start.min(span.end), span.start.max(span.end));
131        return c.sub(a, b, None);
132    }
133    let field = plane_field(curve)?;
134    let n = 128;
135    let at = |i: usize| {
136        evaluate2(
137            curve,
138            span.start + (span.end - span.start) * i as Scalar / n as Scalar,
139        )
140        .ok()
141    };
142    let samples: Vec<Point2> = (0..=n).map(at).collect::<Option<_>>()?;
143    let (mut lo, mut hi) = (samples[0], samples[0]);
144    for p in &samples {
145        lo = lo.min(*p);
146        hi = hi.max(*p);
147    }
148    let pad = (hi - lo) * 0.1 + Vec2::splat(1e-3 * (1.0 + (hi - lo).length()));
149    let curves = crate::implicit_trace::trace(
150        &field,
151        axiolid_curve::implicit::Cell {
152            lo: lo - pad,
153            hi: hi + pad,
154        },
155        crate::implicit_trace::Periodic { u: false, v: false },
156    )
157    .ok()?;
158    let closed = (samples[0] - samples[n]).length() <= 1e-12 * (1.0 + samples[0].length());
159    crate::implicit_ops::extract_stretch(
160        &curves,
161        (false, false),
162        samples[0],
163        [samples[n / 3], samples[2 * n / 3]],
164        samples[n],
165        closed,
166    )
167}
168
169/// Whether `t` lies in `span`, either way round.
170fn within(t: Scalar, span: Interval) -> bool {
171    let (lo, hi) = (span.start.min(span.end), span.start.max(span.end));
172    let slack = 1e-9 * (1.0 + lo.abs().max(hi.abs()));
173    t >= lo - slack && t <= hi + slack
174}
175
176/// Where two plane curves meet, each over its span.
177///
178/// Lines, circles, ellipses, sinusoids, quadratic and angle graphs and
179/// implicit curves, in any pairing. Hits are ordered along the first curve.
180///
181/// # Errors
182///
183/// `UnsupportedCurve` when neither curve has a field (a B-spline or a
184/// lifted pcurve: B-spline pairs take the certified NURBS tier), or when a
185/// piece's trace is refused.
186pub fn section_curve_curve_intersection2(
187    first: &Curve2,
188    first_span: Interval,
189    second: &Curve2,
190    second_span: Interval,
191    tolerance: Tolerance,
192) -> Result<Vec<CurveCurveHit>, ExactCurveRefusal> {
193    // Trace one piece, and find the other's field's roots along it.
194    let (traced, other, swapped) = match (plane_field(first), plane_field(second)) {
195        (_, Some(f)) => (traced_piece(first, first_span), f, false),
196        (Some(f), None) => (traced_piece(second, second_span), f, true),
197        (None, None) => return Err(ExactCurveRefusal::UnsupportedCurve),
198    };
199    let traced = match traced {
200        Some(t) => t,
201        None if plane_field(if swapped { second } else { first }).is_none() => {
202            return Err(ExactCurveRefusal::UnsupportedCurve)
203        }
204        None => return Err(ExactCurveRefusal::UnsupportedCurve),
205    };
206    let roots =
207        crate::implicit_ops::roots_along(&traced, &other, traced.closure(false, false).is_some())
208            .map_err(|_| ExactCurveRefusal::Undecided)?;
209    let mut hits = Vec::new();
210    for (t, multiplicity) in roots {
211        let Some(p) = traced.point(t) else { continue };
212        let (a, b) = match (locate2(first, p, tolerance), locate2(second, p, tolerance)) {
213            (Ok(a), Ok(b)) => (a, b),
214            _ => continue,
215        };
216        if !within(a, first_span) || !within(b, second_span) {
217            continue;
218        }
219        hits.push(CurveCurveHit {
220            first: a,
221            second: b,
222            point: Point3::new(p.x, p.y, 0.0),
223            multiplicity,
224        });
225    }
226    hits.sort_by(|x, y| x.first.total_cmp(&y.first));
227    hits.dedup_by(|x, y| (x.first - y.first).abs() <= 1e-9 * (1.0 + x.first.abs()));
228    Ok(hits)
229}
230
231/// The surface a space curve lies on that best separates it from others:
232/// a traced section's carrier, a conic's plane; a line needs two planes.
233fn carriers_of(curve: &Curve3) -> Option<Vec<Surface>> {
234    let plane = |origin: Point3, z: Vec3| -> Surface {
235        let z = z.normalize();
236        let helper = if z.x.abs() < 0.9 { Vec3::X } else { Vec3::Y };
237        let x = helper.cross(z).normalize();
238        Surface::Plane(Plane {
239            frame: Frame3 {
240                origin,
241                x,
242                y: z.cross(x),
243                z,
244            },
245        })
246    };
247    Some(match curve {
248        Curve3::Line(l) => {
249            let d = l.direction.normalize();
250            let helper = if d.x.abs() < 0.9 { Vec3::X } else { Vec3::Y };
251            let n = d.cross(helper).normalize();
252            vec![plane(l.origin, n), plane(l.origin, d.cross(n))]
253        }
254        Curve3::Circle(c) => vec![plane(c.frame.origin, c.frame.x.cross(c.frame.y))],
255        Curve3::Ellipse(e) => vec![plane(e.frame.origin, e.frame.x.cross(e.frame.y))],
256        Curve3::ImplicitSection(s) => vec![surface_of(&s.carrier)?],
257        Curve3::RuledSection(r) => vec![surface_of(&Carrier::Ruled(r.carrier))?],
258        Curve3::TorusSection(t) => vec![surface_of(&Carrier::Torus(t.torus))?],
259        _ => return None,
260    })
261}
262
263/// A carrier as the surface it copies.
264pub(crate) fn surface_of(carrier: &Carrier) -> Option<Surface> {
265    use axiolid_surface::{Cone, Cylinder, EllipticalCylinder, Sphere, Torus};
266    Some(match carrier {
267        Carrier::Plane(f) => Surface::Plane(Plane { frame: *f }),
268        Carrier::Ruled(k) => {
269            if k.slope != 0.0 {
270                if k.x_radius != k.y_radius {
271                    return None;
272                }
273                Surface::Cone(Cone {
274                    frame: k.frame,
275                    radius: k.x_radius,
276                    semi_angle: k.slope.atan(),
277                })
278            } else if k.x_radius == k.y_radius {
279                Surface::Cylinder(Cylinder {
280                    frame: k.frame,
281                    radius: k.x_radius,
282                })
283            } else {
284                Surface::EllipticalCylinder(EllipticalCylinder {
285                    frame: k.frame,
286                    semi_axis_x: k.x_radius,
287                    semi_axis_y: k.y_radius,
288                })
289            }
290        }
291        Carrier::Sphere { frame, radius } => Surface::Sphere(Sphere {
292            frame: *frame,
293            radius: *radius,
294        }),
295        Carrier::Torus(t) => Surface::Torus(Torus {
296            frame: t.frame,
297            major_radius: t.major_radius,
298            minor_radius: t.minor_radius,
299        }),
300        Carrier::Spline(b) => Surface::BSpline((**b).clone()),
301    })
302}
303
304/// Where two space curves meet, each over its span, within `tolerance`.
305///
306/// The first is intersected with each surface the second lies on
307/// (certified curve/surface intersection), and a point is kept where it
308/// lies on the second curve within `tolerance` and inside both spans.
309///
310/// # Errors
311///
312/// `UnsupportedCurve` for a family neither routine handles.
313pub fn section_curve_curve_intersection3(
314    first: &Curve3,
315    first_span: Interval,
316    second: &Curve3,
317    second_span: Interval,
318    tolerance: Tolerance,
319) -> Result<Vec<CurveCurveHit>, ExactCurveRefusal> {
320    let carriers = carriers_of(second).ok_or(ExactCurveRefusal::UnsupportedCurve)?;
321    let mut hits = Vec::new();
322    for carrier in &carriers {
323        let found = match first {
324            Curve3::Line(_) | Curve3::Circle(_) | Curve3::Ellipse(_) => {
325                crate::exact_curve_intersection::exact_curve_surface_intersection(first, carrier)
326            }
327            _ => {
328                crate::implicit_ops::section_curve_surface_intersection(first, first_span, carrier)
329            }
330        };
331        let points = match found {
332            Ok(ExactCurveIntersection::Points(points)) => points,
333            // The first curve lies on this carrier: the other carriers (a
334            // line's second plane) or the membership test decide.
335            Ok(_) => continue,
336            Err(e) => return Err(e),
337        };
338        for hit in points {
339            let Ok(b) = locate3(second, hit.point, tolerance) else {
340                continue;
341            };
342            let Ok(on) = evaluate3(second, b) else {
343                continue;
344            };
345            if (on - hit.point).length() > tolerance.linear() {
346                continue;
347            }
348            let a = hit.parameter.approx();
349            if !within(a, first_span) || !within(b, second_span) {
350                continue;
351            }
352            hits.push(CurveCurveHit {
353                first: a,
354                second: b,
355                point: hit.point,
356                multiplicity: hit.multiplicity,
357            });
358        }
359    }
360    hits.sort_by(|x, y| x.first.total_cmp(&y.first));
361    hits.dedup_by(|x, y| (x.first - y.first).abs() <= 1e-9 * (1.0 + x.first.abs()));
362    Ok(hits)
363}
axiolid_nurbs/section_curves.rs

axiolid_nurbs/
section_curves.rs
