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

axiolid_nurbs/
section_curves.rs
