axiolid_nurbs/
extrema.rs

1//! Extrema between points, curves and surfaces (#119, B7): the smallest
2//! distance between two pieces of geometry, certified.
3//!
4//! A piece is a point, a curve over a span, or a surface over a box of its
5//! parameters. The search is branch and bound over the pieces' parameter
6//! ranges. Each sub-range has a box certain to hold its image:
7//! - lines and patches of a plane map affinely;
8//! - circles, ellipses, cylinders, cones, spheres and tori have exact
9//!   ranges of their harmonics;
10//! - B-splines are held by the control hull of the restricted Bezier
11//!   pieces (rational weights by interval division);
12//! - traced section curves by their certified cells and their carrier's
13//!   image of them.
14//!
15//! The distance between two boxes bounds the pair from below. The distance
16//! between two evaluated points bounds it from above, and those points
17//! witness it. The search stops when the bounds are within the accuracy
18//! asked for. Nothing is sampled on trust: a pair is discarded only when its
19//! lower bound exceeds the best distance found.
20
21use axiolid_core::{Interval, Point2, Point3, Scalar, Vec3};
22use axiolid_curve::implicit::{bound_simple, partial, Cell};
23use axiolid_curve::{Axis, Carrier, Curve3, ImplicitCurve2};
24use axiolid_evaluate::evaluate3;
25use axiolid_surface::Surface;
26use core::f64::consts::{PI, TAU};
27
28use crate::field::carrier_of;
29
30/// One piece of geometry.
31#[derive(Debug, Clone, Copy)]
32pub enum Piece<'a> {
33    /// A point.
34    Point(Point3),
35    /// A curve over a span of its parameter.
36    Curve {
37        /// The curve.
38        curve: &'a Curve3,
39        /// The span.
40        span: Interval,
41    },
42    /// A surface over a box of its parameters.
43    Patch {
44        /// The surface.
45        surface: &'a Surface,
46        /// Lower corner of the parameter box.
47        lo: Point2,
48        /// Upper corner.
49        hi: Point2,
50    },
51}
52
53/// The smallest distance between two pieces, bracketed, with points on
54/// each at the upper bound.
55#[derive(Debug, Clone, Copy, PartialEq)]
56pub struct Extremum {
57    /// Certain lower bound on the distance.
58    pub lower: Scalar,
59    /// Distance between the witnesses: an upper bound.
60    pub upper: Scalar,
61    /// Point on the first piece.
62    pub on_first: Point3,
63    /// Point on the second piece.
64    pub on_second: Point3,
65    /// Its parameters on the first piece (`t`, or `(u, v)`; `NaN` for a
66    /// point).
67    pub first: Point2,
68    /// Its parameters on the second piece.
69    pub second: Point2,
70}
71
72/// Why an extremum could not be found.
73#[derive(Debug, Clone, Copy, PartialEq, Eq)]
74#[non_exhaustive]
75pub enum ExtremaRefusal {
76    /// A curve or surface family this module has no image bounds for.
77    Unsupported,
78    /// A piece could not be evaluated.
79    Evaluation,
80    /// The search exceeded its work budget before reaching the accuracy.
81    Budget,
82}
83
84/// An axis-aligned box.
85#[derive(Debug, Clone, Copy)]
86struct Aabb {
87    lo: Vec3,
88    hi: Vec3,
89}
90
91impl Aabb {
92    fn point(p: Point3) -> Self {
93        Self { lo: p, hi: p }
94    }
95
96    fn join(self, o: Self) -> Self {
97        Self {
98            lo: self.lo.min(o.lo),
99            hi: self.hi.max(o.hi),
100        }
101    }
102
103    fn distance(&self, o: &Self) -> Scalar {
104        let gap = (self.lo - o.hi).max(o.lo - self.hi).max(Vec3::ZERO);
105        gap.length()
106    }
107
108    fn diameter(&self) -> Scalar {
109        (self.hi - self.lo).length()
110    }
111
112    fn widen(self, by: Scalar) -> Self {
113        Self {
114            lo: self.lo - Vec3::splat(by),
115            hi: self.hi + Vec3::splat(by),
116        }
117    }
118}
119
120/// An interval of reals, with the arithmetic the bounds need.
121#[derive(Debug, Clone, Copy)]
122struct Iv(Scalar, Scalar);
123
124impl Iv {
125    fn add(self, o: Iv) -> Iv {
126        Iv(self.0 + o.0, self.1 + o.1)
127    }
128
129    fn mul(self, o: Iv) -> Iv {
130        let p = [self.0 * o.0, self.0 * o.1, self.1 * o.0, self.1 * o.1];
131        Iv(
132            p.iter().copied().fold(Scalar::INFINITY, Scalar::min),
133            p.iter().copied().fold(Scalar::NEG_INFINITY, Scalar::max),
134        )
135    }
136
137    fn scale(self, c: Scalar) -> Iv {
138        let (a, b) = (self.0 * c, self.1 * c);
139        Iv(a.min(b), a.max(b))
140    }
141}
142
143fn cos_iv(a: Scalar, b: Scalar) -> Iv {
144    if b - a >= TAU {
145        return Iv(-1.0, 1.0);
146    }
147    let (ca, cb) = (a.cos(), b.cos());
148    let (mut lo, mut hi) = (ca.min(cb), ca.max(cb));
149    if ((a / TAU).ceil()) * TAU <= b {
150        hi = 1.0;
151    }
152    if (((a - PI) / TAU).ceil()) * TAU + PI <= b {
153        lo = -1.0;
154    }
155    Iv(lo, hi)
156}
157
158fn sin_iv(a: Scalar, b: Scalar) -> Iv {
159    cos_iv(a - 0.5 * PI, b - 0.5 * PI)
160}
161
162/// `origin + sum_k axis_k * factor_k` over intervals of the factors.
163fn combine(origin: Point3, terms: &[(Vec3, Iv)]) -> Aabb {
164    let mut lo = origin;
165    let mut hi = origin;
166    for (axis, factor) in terms {
167        for k in 0..3 {
168            let c = factor.scale(axis[k]);
169            lo[k] += c.0;
170            hi[k] += c.1;
171        }
172    }
173    Aabb { lo, hi }
174}
175
176/// A box holding the carrier's points over a parameter box.
177fn carrier_box(carrier: &Carrier, lo: Point2, hi: Point2) -> Option<Aabb> {
178    let (cu, su) = (cos_iv(lo.x, hi.x), sin_iv(lo.x, hi.x));
179    let bx = match carrier {
180        Carrier::Plane(f) => {
181            let mut b = Aabb::point(f.origin + f.x * lo.x + f.y * lo.y);
182            for (u, v) in [(hi.x, lo.y), (lo.x, hi.y), (hi.x, hi.y)] {
183                b = b.join(Aabb::point(f.origin + f.x * u + f.y * v));
184            }
185            b
186        }
187        Carrier::Ruled(k) => {
188            let v = Iv(lo.y, hi.y);
189            let rx = v.scale(k.slope).add(Iv(k.x_radius, k.x_radius));
190            let ry = v.scale(k.slope).add(Iv(k.y_radius, k.y_radius));
191            combine(
192                k.frame.origin,
193                &[
194                    (k.frame.x, rx.mul(cu)),
195                    (k.frame.y, ry.mul(su)),
196                    (k.frame.z, v),
197                ],
198            )
199        }
200        Carrier::Sphere { frame, radius } => {
201            let (cv, sv) = (cos_iv(lo.y, hi.y), sin_iv(lo.y, hi.y));
202            let r = Iv(*radius, *radius);
203            combine(
204                frame.origin,
205                &[
206                    (frame.x, r.mul(cv).mul(cu)),
207                    (frame.y, r.mul(cv).mul(su)),
208                    (frame.z, r.mul(sv)),
209                ],
210            )
211        }
212        Carrier::Torus(t) => {
213            let (cv, sv) = (cos_iv(lo.y, hi.y), sin_iv(lo.y, hi.y));
214            let ring = cv
215                .scale(t.minor_radius)
216                .add(Iv(t.major_radius, t.major_radius));
217            combine(
218                t.frame.origin,
219                &[
220                    (t.frame.x, ring.mul(cu)),
221                    (t.frame.y, ring.mul(su)),
222                    (t.frame.z, sv.scale(t.minor_radius)),
223                ],
224            )
225        }
226        Carrier::Spline(b) => return spline_box(b, lo, hi),
227    };
228    // Rounding in the sums.
229    let size = bx.lo.abs().max(bx.hi.abs()).max_element();
230    Some(bx.widen(32.0 * Scalar::EPSILON * (1.0 + size)))
231}
232
233/// A box holding a B-spline surface over a parameter box: each coordinate's
234/// homogeneous field bounded by its restricted Bernstein coefficients,
235/// divided by the weight's bounds.
236#[allow(clippy::needless_range_loop)]
237fn spline_box(b: &axiolid_curve::BSplineSurface, lo: Point2, hi: Point2) -> Option<Aabb> {
238    let fields = crate::spline_field::homogeneous_fields(b)?;
239    let cell = Cell { lo, hi };
240    let w = bound_simple(&fields[3], &cell);
241    if w.lo <= 0.0 {
242        return None;
243    }
244    let mut out = Aabb {
245        lo: Vec3::ZERO,
246        hi: Vec3::ZERO,
247    };
248    for k in 0..3 {
249        let x = bound_simple(&fields[k], &cell);
250        let q = [x.lo / w.lo, x.lo / w.hi, x.hi / w.lo, x.hi / w.hi];
251        out.lo[k] = q.iter().copied().fold(Scalar::INFINITY, Scalar::min);
252        out.hi[k] = q.iter().copied().fold(Scalar::NEG_INFINITY, Scalar::max);
253    }
254    Some(out)
255}
256
257/// The box `[free0, free1] x [...]` of carrier parameters certain to hold
258/// an implicit curve's cell for `s` in `[s0, s1]`.
259fn cell_hull(
260    curve: &ImplicitCurve2,
261    index: usize,
262    s0: Scalar,
263    s1: Scalar,
264) -> Option<(Point2, Point2)> {
265    let cell = &curve.cells[index];
266    let free = |s: Scalar| cell.from + (cell.to - cell.from) * s;
267    let (f0, f1) = (free(s0), free(s1));
268    let (f_lo, f_hi) = (f0.min(f1), f0.max(f1));
269    let place = |fr: Scalar, w: Scalar| match cell.axis {
270        Axis::U => Point2::new(fr, w),
271        Axis::V => Point2::new(w, fr),
272    };
273    let w0 = match cell.axis {
274        Axis::U => curve.point(index as Scalar + s0)?.y,
275        Axis::V => curve.point(index as Scalar + s0)?.x,
276    };
277    // A bridge is straight: the box of its ends.
278    // A bridge lies in its bracket about the straight line.
279    if cell.bridge.is_some() {
280        let (w_lo, w_hi) = cell.solved_range();
281        let (a, b) = (place(f_lo, w_lo), place(f_hi, w_hi));
282        return Some((a.min(b), a.max(b)));
283    }
284    let (d_free, d_solved) = match cell.axis {
285        Axis::U => (partial(&curve.field, true), partial(&curve.field, false)),
286        Axis::V => (partial(&curve.field, false), partial(&curve.field, true)),
287    };
288    let whole = {
289        let (a, b) = (place(f_lo, cell.low), place(f_hi, cell.high));
290        Cell {
291            lo: a.min(b),
292            hi: a.max(b),
293        }
294    };
295    let fr = bound_simple(&d_free, &whole);
296    let so = bound_simple(&d_solved, &whole);
297    let (w_lo, w_hi) = if so.straddles_zero() {
298        (cell.low, cell.high)
299    } else {
300        let floor = so.lo.abs().min(so.hi.abs());
301        let reach = fr.lo.abs().max(fr.hi.abs()) / floor * (f_hi - f_lo);
302        ((w0 - reach).max(cell.low), (w0 + reach).min(cell.high))
303    };
304    let (a, b) = (place(f_lo, w_lo), place(f_hi, w_hi));
305    Some((a.min(b), a.max(b)))
306}
307
308/// The range of `a cos t + b sin t` over `[t0, t1]`, exactly.
309fn harmonic(a: Scalar, b: Scalar, t0: Scalar, t1: Scalar) -> Iv {
310    let r = a.hypot(b);
311    if r == 0.0 {
312        return Iv(0.0, 0.0);
313    }
314    let phi = b.atan2(a);
315    cos_iv(t0 - phi, t1 - phi).scale(r)
316}
317
318/// The range of `n . P` over the carrier's points on a parameter box: the
319/// direct bound intersected with the mean-value form
320/// `f(c) + [f_u] [-h_u, h_u] + [f_v] [-h_v, h_v]`, which is second order
321/// where `n` is normal to the surface -- exactly where a closest point is,
322/// so the search there needs boxes only as small as the accuracy's square
323/// root.
324fn carrier_projection(carrier: &Carrier, n: Vec3, lo: Point2, hi: Point2) -> Option<Iv> {
325    let direct = carrier_projection_direct(carrier, n, lo, hi)?;
326    let Some((fu, fv)) = carrier_gradient(carrier, n, lo, hi) else {
327        return Some(direct);
328    };
329    let c = (lo + hi) * 0.5;
330    let (hu, hv) = (0.5 * (hi.x - lo.x), 0.5 * (hi.y - lo.y));
331    let jet = carrier.jet(c.x, c.y);
332    let f = n.dot(jet.point);
333    let m = fu.scale(1.0).mul(Iv(-hu, hu)).add(fv.mul(Iv(-hv, hv)));
334    let slack = 32.0 * Scalar::EPSILON * (1.0 + jet.point.length());
335    let mean = Iv(f + m.0 - slack, f + m.1 + slack);
336    Some(Iv(direct.0.max(mean.0), direct.1.min(mean.1)))
337}
338
339/// Ranges of `n . P_u` and `n . P_v` over a parameter box.
340fn carrier_gradient(carrier: &Carrier, n: Vec3, lo: Point2, hi: Point2) -> Option<(Iv, Iv)> {
341    let (cu, su) = (cos_iv(lo.x, hi.x), sin_iv(lo.x, hi.x));
342    Some(match carrier {
343        Carrier::Plane(f) => {
344            let (a, b) = (n.dot(f.x), n.dot(f.y));
345            (Iv(a, a), Iv(b, b))
346        }
347        Carrier::Ruled(k) => {
348            let f = &k.frame;
349            let (a, b, c) = (n.dot(f.x), n.dot(f.y), n.dot(f.z));
350            let v = Iv(lo.y, hi.y);
351            let rx = v.scale(k.slope).add(Iv(k.x_radius, k.x_radius));
352            let ry = v.scale(k.slope).add(Iv(k.y_radius, k.y_radius));
353            // P_u = -rx sin u X + ry cos u Y; P_v = s cos u X + s sin u Y + Z.
354            let pu = rx.mul(su.scale(-a)).add(ry.mul(cu.scale(b)));
355            let pv = harmonic(a * k.slope, b * k.slope, lo.x, hi.x).add(Iv(c, c));
356            (pu, pv)
357        }
358        Carrier::Sphere { frame, radius } => {
359            let (a, b, c) = (n.dot(frame.x), n.dot(frame.y), n.dot(frame.z));
360            let (cv, sv) = (cos_iv(lo.y, hi.y), sin_iv(lo.y, hi.y));
361            // P_u = r cos v (-sin u X + cos u Y); P_v = r (-sin v (cos u X +
362            // sin u Y) + cos v Z).
363            let pu = cv.mul(harmonic(b, -a, lo.x, hi.x)).scale(*radius);
364            let pv = sv
365                .mul(harmonic(a, b, lo.x, hi.x))
366                .scale(-1.0)
367                .add(cv.scale(c))
368                .scale(*radius);
369            (pu, pv)
370        }
371        Carrier::Torus(t) => {
372            let f = &t.frame;
373            let (a, b, c) = (n.dot(f.x), n.dot(f.y), n.dot(f.z));
374            let (cv, sv) = (cos_iv(lo.y, hi.y), sin_iv(lo.y, hi.y));
375            let r = t.minor_radius;
376            let ring = cv.scale(r).add(Iv(t.major_radius, t.major_radius));
377            let pu = ring.mul(harmonic(b, -a, lo.x, hi.x));
378            let pv = sv
379                .mul(harmonic(a, b, lo.x, hi.x))
380                .scale(-r)
381                .add(cv.scale(c * r));
382            (pu, pv)
383        }
384        Carrier::Spline(_) => return None,
385    })
386}
387
388/// The direct bound of [`carrier_projection`].
389fn carrier_projection_direct(carrier: &Carrier, n: Vec3, lo: Point2, hi: Point2) -> Option<Iv> {
390    let pad = |iv: Iv, scale: Scalar| {
391        let m = 32.0 * Scalar::EPSILON * (1.0 + scale);
392        Iv(iv.0 - m, iv.1 + m)
393    };
394    Some(match carrier {
395        Carrier::Plane(f) => {
396            let o = n.dot(f.origin);
397            let (a, b) = (n.dot(f.x), n.dot(f.y));
398            let (u, v) = (Iv(lo.x, hi.x).scale(a), Iv(lo.y, hi.y).scale(b));
399            pad(Iv(o, o).add(u).add(v), o.abs() + a.abs() + b.abs())
400        }
401        Carrier::Ruled(k) => {
402            let f = &k.frame;
403            let (a, b, c) = (n.dot(f.x), n.dot(f.y), n.dot(f.z));
404            let o = n.dot(f.origin);
405            let v = Iv(lo.y, hi.y);
406            let iv = if k.slope == 0.0 && k.x_radius == k.y_radius {
407                // A cylinder: separable, exact.
408                harmonic(a * k.x_radius, b * k.y_radius, lo.x, hi.x).add(v.scale(c))
409            } else {
410                let rx = v.scale(k.slope).add(Iv(k.x_radius, k.x_radius));
411                let ry = v.scale(k.slope).add(Iv(k.y_radius, k.y_radius));
412                rx.mul(cos_iv(lo.x, hi.x).scale(a))
413                    .add(ry.mul(sin_iv(lo.x, hi.x).scale(b)))
414                    .add(v.scale(c))
415            };
416            pad(
417                Iv(o, o).add(iv),
418                o.abs() + k.x_radius + k.y_radius + v.1.abs(),
419            )
420        }
421        Carrier::Sphere { frame, radius } => {
422            let (a, b, c) = (n.dot(frame.x), n.dot(frame.y), n.dot(frame.z));
423            let o = n.dot(frame.origin);
424            let h = harmonic(a, b, lo.x, hi.x);
425            let iv = cos_iv(lo.y, hi.y)
426                .mul(h)
427                .add(sin_iv(lo.y, hi.y).scale(c))
428                .scale(*radius);
429            pad(Iv(o, o).add(iv), o.abs() + radius)
430        }
431        Carrier::Torus(t) => {
432            let f = &t.frame;
433            let (a, b, c) = (n.dot(f.x), n.dot(f.y), n.dot(f.z));
434            let o = n.dot(f.origin);
435            let h = harmonic(a, b, lo.x, hi.x);
436            let ring = cos_iv(lo.y, hi.y)
437                .scale(t.minor_radius)
438                .add(Iv(t.major_radius, t.major_radius));
439            let iv = ring
440                .mul(h)
441                .add(sin_iv(lo.y, hi.y).scale(c * t.minor_radius));
442            pad(Iv(o, o).add(iv), o.abs() + t.major_radius + t.minor_radius)
443        }
444        Carrier::Spline(b) => {
445            let bx = spline_box(b, lo, hi)?;
446            project_box(&bx, n)
447        }
448    })
449}
450
451/// The range of `n . P` over a box.
452fn project_box(b: &Aabb, n: Vec3) -> Iv {
453    let mut lo = 0.0;
454    let mut hi = 0.0;
455    for k in 0..3 {
456        let (x, y) = (n[k] * b.lo[k], n[k] * b.hi[k]);
457        lo += x.min(y);
458        hi += x.max(y);
459    }
460    Iv(lo, hi)
461}
462
463/// The range of `n . C(t)` over `[t0, t1]` of a curve.
464fn curve_projection(curve: &Curve3, n: Vec3, t0: Scalar, t1: Scalar) -> Option<Iv> {
465    let (a, b) = (t0.min(t1), t0.max(t1));
466    match curve {
467        Curve3::Line(l) => {
468            let (p, q) = (
469                n.dot(l.origin + l.direction * a),
470                n.dot(l.origin + l.direction * b),
471            );
472            Some(Iv(p.min(q), p.max(q)))
473        }
474        // Exact: one harmonic.
475        Curve3::Circle(c) => {
476            let o = n.dot(c.frame.origin);
477            let h = harmonic(
478                n.dot(c.frame.x) * c.radius,
479                n.dot(c.frame.y) * c.radius,
480                a,
481                b,
482            );
483            Some(Iv(o, o).add(h))
484        }
485        Curve3::Ellipse(e) => {
486            let o = n.dot(e.frame.origin);
487            let h = harmonic(
488                n.dot(e.frame.x) * e.semi_axis_x,
489                n.dot(e.frame.y) * e.semi_axis_y,
490                a,
491                b,
492            );
493            Some(Iv(o, o).add(h))
494        }
495        Curve3::ImplicitSection(s) => {
496            let mut out: Option<Iv> = None;
497            for index in 0..s.curve.cells.len() {
498                let (c0, c1) = (index as Scalar, index as Scalar + 1.0);
499                let (lo, hi) = (a.max(c0), b.min(c1));
500                if hi < lo {
501                    continue;
502                }
503                let (p, q) = cell_hull(&s.curve, index, lo - c0, hi - c0)?;
504                let iv = carrier_projection(&s.carrier, n, p, q)?;
505                out = Some(out.map_or(iv, |o: Iv| Iv(o.0.min(iv.0), o.1.max(iv.1))));
506            }
507            out
508        }
509        _ => curve_box(curve, a, b).map(|bx| project_box(&bx, n)),
510    }
511}
512
513/// A box holding a curve's image over `[t0, t1]`.
514fn curve_box(curve: &Curve3, t0: Scalar, t1: Scalar) -> Option<Aabb> {
515    let (a, b) = (t0.min(t1), t0.max(t1));
516    match curve {
517        Curve3::Line(l) => Some(
518            Aabb::point(l.origin + l.direction * a).join(Aabb::point(l.origin + l.direction * b)),
519        ),
520        Curve3::Circle(c) => Some(combine(
521            c.frame.origin,
522            &[
523                (c.frame.x, cos_iv(a, b).scale(c.radius)),
524                (c.frame.y, sin_iv(a, b).scale(c.radius)),
525            ],
526        )),
527        Curve3::Ellipse(e) => Some(combine(
528            e.frame.origin,
529            &[
530                (e.frame.x, cos_iv(a, b).scale(e.semi_axis_x)),
531                (e.frame.y, sin_iv(a, b).scale(e.semi_axis_y)),
532            ],
533        )),
534        Curve3::ImplicitSection(s) => {
535            let mut out: Option<Aabb> = None;
536            for index in 0..s.curve.cells.len() {
537                let (c0, c1) = (index as Scalar, index as Scalar + 1.0);
538                let (lo, hi) = (a.max(c0), b.min(c1));
539                if hi < lo {
540                    continue;
541                }
542                let (p, q) = cell_hull(&s.curve, index, lo - c0, hi - c0)?;
543                let bx = carrier_box(&s.carrier, p, q)?;
544                out = Some(out.map_or(bx, |o| o.join(bx)));
545            }
546            out
547        }
548        Curve3::BSpline(_) => {
549            // Its control hull over the span, by the curve's own Bezier
550            // pieces: bounded conservatively by the whole control polygon
551            // restricted to the span's knot range.
552            crate::spline_field::curve_hull(curve, a, b).map(|(lo, hi)| Aabb { lo, hi })
553        }
554        _ => None,
555    }
556}
557
558/// A piece reduced to what the search needs: a parameter range, how to
559/// bound its image and where its points are.
560#[derive(Debug, Clone)]
561enum Part {
562    Point(Point3),
563    Curve {
564        curve: Curve3,
565        lo: Scalar,
566        hi: Scalar,
567    },
568    Patch {
569        surface: Surface,
570        carrier: Carrier,
571        lo: Point2,
572        hi: Point2,
573    },
574}
575
576impl Part {
577    fn of(piece: &Piece<'_>) -> Result<Part, ExtremaRefusal> {
578        Ok(match piece {
579            Piece::Point(p) => Part::Point(*p),
580            Piece::Curve { curve, span } => {
581                // The ADR 0076 graphs are bounded through their traced form.
582                let curve = match curve {
583                    Curve3::RuledSection(_) | Curve3::TorusSection(_) => Curve3::ImplicitSection(
584                        crate::implicit_ops::implicit_view(curve, *span)
585                            .ok_or(ExtremaRefusal::Unsupported)?,
586                    ),
587                    other => (*other).clone(),
588                };
589                let (lo, hi) = match (&curve, piece) {
590                    (
591                        Curve3::ImplicitSection(s),
592                        Piece::Curve {
593                            curve: original, ..
594                        },
595                    ) if !matches!(original, Curve3::ImplicitSection(_)) => (0.0, s.curve.end()),
596                    _ => (span.start.min(span.end), span.start.max(span.end)),
597                };
598                Part::Curve { curve, lo, hi }
599            }
600            Piece::Patch { surface, lo, hi } => Part::Patch {
601                surface: (*surface).clone(),
602                carrier: carrier_of(surface).ok_or(ExtremaRefusal::Unsupported)?,
603                lo: *lo,
604                hi: *hi,
605            },
606        })
607    }
608}
609
610/// A sub-range of a part.
611#[derive(Debug, Clone, Copy)]
612enum Range2 {
613    Point,
614    Curve(Scalar, Scalar),
615    Patch(Point2, Point2),
616}
617
618impl Range2 {
619    fn whole(part: &Part) -> Range2 {
620        match part {
621            Part::Point(_) => Range2::Point,
622            Part::Curve { lo, hi, .. } => Range2::Curve(*lo, *hi),
623            Part::Patch { lo, hi, .. } => Range2::Patch(*lo, *hi),
624        }
625    }
626
627    fn bounds(&self, part: &Part) -> Option<Aabb> {
628        match (self, part) {
629            (Range2::Point, Part::Point(p)) => Some(Aabb::point(*p)),
630            (Range2::Curve(a, b), Part::Curve { curve, .. }) => curve_box(curve, *a, *b),
631            (Range2::Patch(lo, hi), Part::Patch { carrier, .. }) => carrier_box(carrier, *lo, *hi),
632            _ => None,
633        }
634    }
635
636    /// The range of `n . P` over the part's points in this range.
637    fn projection(&self, part: &Part, n: Vec3) -> Option<Iv> {
638        match (self, part) {
639            (Range2::Point, Part::Point(p)) => {
640                let x = n.dot(*p);
641                Some(Iv(x, x))
642            }
643            (Range2::Curve(a, b), Part::Curve { curve, .. }) => curve_projection(curve, n, *a, *b),
644            (Range2::Patch(lo, hi), Part::Patch { carrier, .. }) => {
645                carrier_projection(carrier, n, *lo, *hi)
646            }
647            _ => None,
648        }
649    }
650
651    /// The part's point in the middle of the range, and its parameters.
652    fn sample(&self, part: &Part) -> Option<(Point3, Point2)> {
653        match (self, part) {
654            (Range2::Point, Part::Point(p)) => Some((*p, Point2::splat(Scalar::NAN))),
655            (Range2::Curve(a, b), Part::Curve { curve, .. }) => {
656                let t = 0.5 * (a + b);
657                Some((evaluate3(curve, t).ok()?, Point2::new(t, Scalar::NAN)))
658            }
659            (Range2::Patch(lo, hi), Part::Patch { surface, .. }) => {
660                let c = (*lo + *hi) * 0.5;
661                Some((
662                    axiolid_evaluate::surface::evaluate(surface, c.x, c.y).ok()?,
663                    c,
664                ))
665            }
666            _ => None,
667        }
668    }
669
670    fn split(&self) -> Vec<Range2> {
671        match self {
672            Range2::Point => vec![Range2::Point],
673            Range2::Curve(a, b) => {
674                let m = 0.5 * (a + b);
675                vec![Range2::Curve(*a, m), Range2::Curve(m, *b)]
676            }
677            Range2::Patch(lo, hi) => {
678                let d = *hi - *lo;
679                if d.x >= d.y {
680                    let m = 0.5 * (lo.x + hi.x);
681                    vec![
682                        Range2::Patch(*lo, Point2::new(m, hi.y)),
683                        Range2::Patch(Point2::new(m, lo.y), *hi),
684                    ]
685                } else {
686                    let m = 0.5 * (lo.y + hi.y);
687                    vec![
688                        Range2::Patch(*lo, Point2::new(hi.x, m)),
689                        Range2::Patch(Point2::new(lo.x, m), *hi),
690                    ]
691                }
692            }
693        }
694    }
695
696    fn is_point(&self) -> bool {
697        matches!(self, Range2::Point)
698    }
699}
700
701/// The smallest distance between two pieces, to within `accuracy`.
702///
703/// # Errors
704///
705/// A family without image bounds, an evaluation failure, or a search that
706/// needs more than its budget of pairs.
707pub fn minimum_distance(
708    first: &Piece<'_>,
709    second: &Piece<'_>,
710    accuracy: Scalar,
711) -> Result<Extremum, ExtremaRefusal> {
712    let (pa, pb) = (Part::of(first)?, Part::of(second)?);
713    let (ra, rb) = (Range2::whole(&pa), Range2::whole(&pb));
714    let sample = |r: &Range2, p: &Part| r.sample(p).ok_or(ExtremaRefusal::Evaluation);
715    let (sa, sb) = (sample(&ra, &pa)?, sample(&rb, &pb)?);
716    let mut best = Extremum {
717        lower: 0.0,
718        upper: (sa.0 - sb.0).length(),
719        on_first: sa.0,
720        on_second: sb.0,
721        first: sa.1,
722        second: sb.1,
723    };
724    let bounds = |r: &Range2, p: &Part| r.bounds(p).ok_or(ExtremaRefusal::Unsupported);
725    // A lower bound for a pair: the gap between their boxes, or better the
726    // gap between their projections on the line joining their boxes'
727    // centres, which closes on the true distance as the pieces shrink.
728    let lower_bound = |x: &Range2, y: &Range2, bx: &Aabb, by: &Aabb| -> Scalar {
729        let boxes = bx.distance(by);
730        let d = (bx.lo + bx.hi) * 0.5 - (by.lo + by.hi) * 0.5;
731        let length = d.length();
732        if length == 0.0 {
733            return boxes;
734        }
735        let n = d / length;
736        match (x.projection(&pa, n), y.projection(&pb, n)) {
737            (Some(px), Some(py)) => boxes.max(px.0 - py.1),
738            _ => boxes,
739        }
740    };
741    // Pending pairs with their lower bounds, the smallest searched first.
742    let (ba, bb) = (bounds(&ra, &pa)?, bounds(&rb, &pb)?);
743    let mut pending: Vec<(Scalar, Range2, Range2)> =
744        vec![(lower_bound(&ra, &rb, &ba, &bb), ra, rb)];
745    let mut work = 0usize;
746    loop {
747        // The smallest lower bound still open.
748        let Some(index) = pending
749            .iter()
750            .enumerate()
751            .min_by(|x, y| x.1 .0.total_cmp(&y.1 .0))
752            .map(|(i, _)| i)
753        else {
754            best.lower = best.upper;
755            return Ok(best);
756        };
757        let (lower, a, b) = pending.swap_remove(index);
758        best.lower = lower.min(best.upper);
759        if best.upper - lower <= accuracy {
760            return Ok(best);
761        }
762        work += 1;
763        if work > 200_000 {
764            return Err(ExtremaRefusal::Budget);
765        }
766        // Split the side whose image is larger.
767        let (da, db) = (bounds(&a, &pa)?.diameter(), bounds(&b, &pb)?.diameter());
768        let (split_first, parts) = if (da >= db && !a.is_point()) || b.is_point() {
769            (true, a.split())
770        } else {
771            (false, b.split())
772        };
773        for part in parts {
774            let (x, y) = if split_first { (part, b) } else { (a, part) };
775            let (bx, by) = (bounds(&x, &pa)?, bounds(&y, &pb)?);
776            let low = lower_bound(&x, &y, &bx, &by);
777            if low > best.upper {
778                continue;
779            }
780            let (p, q) = (sample(&x, &pa)?, sample(&y, &pb)?);
781            let d = (p.0 - q.0).length();
782            if d < best.upper {
783                best.upper = d;
784                best.on_first = p.0;
785                best.on_second = q.0;
786                best.first = p.1;
787                best.second = q.1;
788            }
789            if low <= best.upper {
790                pending.push((low, x, y));
791            }
792        }
793        pending.retain(|(low, _, _)| *low <= best.upper);
794    }
795}
axiolid_nurbs/extrema.rs

axiolid_nurbs/
extrema.rs
