axiolid_inspect/
overlap.rs

1//! Certified volume shared by two closed meshes, and the volume of their
2//! difference (#183).
3//!
4//! # Prisms, not a boolean
5//!
6//! No mesh is built. For a closed, consistently oriented, non-self-
7//! intersecting triangle mesh `A` lying above a plane `z = z0`, a point is
8//! inside `A` exactly when the faces above it, counted `+1` for faces
9//! facing up and `-1` for faces facing down, sum to one. So, almost
10//! everywhere,
11//!
12//! ```text
13//! 1_A = sum over faces f of  s_f * 1{ (x, y) in f', z0 < z < h_f(x, y) }
14//! ```
15//!
16//! with `f'` the face's shadow on the plane, `h_f` the height of its plane
17//! and `s_f` the sign of the shadow's orientation (vertical faces cast no
18//! shadow and drop out). Multiplying two such sums and integrating,
19//!
20//! ```text
21//! vol(A and B) = sum over f in A, g in B of  s_f s_g * integral over f' and g' of (min(h_f, h_g) - z0)
22//! ```
23//!
24//! Each term is a convex polygon -- the overlap of two triangles --
25//! split by the line where the two planes cross, and a linear function
26//! integrated over each piece.
27//!
28//! # Exact decisions, enclosed arithmetic
29//!
30//! Every vertex of every piece before the split is an input vertex or the
31//! crossing of two input edge lines, so which side of a line or of the
32//! plane crossing it lies on is the sign of a polynomial in the input
33//! coordinates: decided by interval arithmetic and else exactly in
34//! dyadics. The numbers -- crossing points, heights, areas -- are then
35//! computed in outward-rounded intervals, so the sum is an interval that
36//! contains the true volume. Touching bodies and shared faces cancel
37//! exactly in the decisions; their volume comes out as a small interval
38//! around zero, and is clamped to be no less than zero.
39//!
40//! # A priori bound per term (#313)
41//!
42//! A face that is vertical but for rounding -- corners that differ by
43//! about `1e-15` across, as exported building models have -- still casts a
44//! shadow, a sliver. Its plane's slope is then near `1e30`, so its height
45//! at a crossing point enclosed to `1e-15` spans `±1e15`, and the sliver's
46//! term, though its true value is tiny, comes out wider than either solid.
47//! So each term is also intersected with a bound that does not depend on
48//! how well its arithmetic is conditioned. Over its own shadow a face's
49//! plane is the triangle itself, so its height lies between its lowest and
50//! highest corner; the term's region lies inside both shadows; and `z0` is
51//! no higher than any corner of either mesh. Hence
52//!
53//! ```text
54//! 0 <= term <= min(area f', area g') * (min(zmax_f, zmax_g) - z0)
55//! ```
56//!
57//! computed outward-rounded. A sliver's term is then no wider than the
58//! sliver's area times its height, and the bound only ever narrows an
59//! interval that already holds the true value, so exactly coplanar faces
60//! cancel as before.
61
62use axiolid_core::{Aabb, Point2, Point3, Tolerance};
63use axiolid_exact::{certify, Arith, SignExpr};
64use axiolid_guarantees::Sign;
65use axiolid_heal::self_intersections;
66use axiolid_mesh::{audit_mesh, TriangleMeshView};
67use axiolid_spatial::{Bvh, SpatialIndex, SpatialItem};
68use core::ops::ControlFlow;
69
70/// A closed interval of volumes.
71#[derive(Debug, Clone, Copy, PartialEq)]
72pub struct VolumeInterval {
73    /// No greater than the true volume, and never negative.
74    pub lower: f64,
75    /// No less than the true volume.
76    pub upper: f64,
77}
78
79impl VolumeInterval {
80    /// Whether `volume` lies in the interval.
81    #[must_use]
82    pub fn contains(&self, volume: f64) -> bool {
83        self.lower <= volume && volume <= self.upper
84    }
85
86    /// `upper - lower`.
87    #[must_use]
88    pub fn width(&self) -> f64 {
89        self.upper - self.lower
90    }
91}
92
93/// Which operand an error is about.
94#[derive(Debug, Clone, Copy, PartialEq, Eq)]
95#[non_exhaustive]
96pub enum Operand {
97    /// The first mesh given.
98    First,
99    /// The second mesh given.
100    Second,
101}
102
103/// Why no volume was bracketed.
104#[derive(Debug, Clone, Copy, PartialEq, Eq)]
105#[non_exhaustive]
106pub enum OverlapError {
107    /// A coordinate is not finite.
108    NonFinite {
109        /// Which mesh.
110        operand: Operand,
111    },
112    /// The mesh is not a closed, consistently wound two-manifold without
113    /// degenerate triangles, so it encloses no definite solid.
114    NotClosed {
115        /// Which mesh.
116        operand: Operand,
117    },
118    /// Two triangles of the mesh intersect other than along shared
119    /// vertices, so points inside it are not well defined.
120    SelfIntersecting {
121        /// Which mesh.
122        operand: Operand,
123        /// The first pair found, lower index first.
124        triangles: [u32; 2],
125    },
126    /// The mesh's enclosed volume is not certified nonzero.
127    NoVolume {
128        /// Which mesh.
129        operand: Operand,
130    },
131}
132
133/// The volume enclosed by a closed mesh, certified. Either winding is
134/// accepted.
135///
136/// # Errors
137///
138/// [`OverlapError`] for a non-finite, open or self-intersecting mesh, or
139/// one whose volume is not certified nonzero.
140pub fn enclosed_volume<M: TriangleMeshView + ?Sized>(
141    mesh: &M,
142) -> Result<VolumeInterval, OverlapError> {
143    let solid = Solid::new(mesh, Operand::First, None)?;
144    Ok(clamp(solid.volume, 0.0, f64::INFINITY))
145}
146
147/// The volume two closed meshes share, certified.
148///
149/// # Errors
150///
151/// [`OverlapError`] naming the operand that is non-finite, open, self-
152/// intersecting or without volume.
153pub fn intersection_volume<A, B>(first: &A, second: &B) -> Result<VolumeInterval, OverlapError>
154where
155    A: TriangleMeshView + ?Sized,
156    B: TriangleMeshView + ?Sized,
157{
158    let (a, b) = solids(first, second)?;
159    Ok(shared(&a, &b))
160}
161
162/// The volume of the first mesh outside the second, certified.
163///
164/// # Errors
165///
166/// As [`intersection_volume`].
167pub fn difference_volume<A, B>(first: &A, second: &B) -> Result<VolumeInterval, OverlapError>
168where
169    A: TriangleMeshView + ?Sized,
170    B: TriangleMeshView + ?Sized,
171{
172    let (a, b) = solids(first, second)?;
173    let both = shared(&a, &b);
174    let rest = a.volume.sub(Iv::new(both.lower, both.upper));
175    Ok(clamp(rest, 0.0, a.volume.hi))
176}
177
178fn solids<A, B>(first: &A, second: &B) -> Result<(Solid, Solid), OverlapError>
179where
180    A: TriangleMeshView + ?Sized,
181    B: TriangleMeshView + ?Sized,
182{
183    let base = lowest(first).min(lowest(second));
184    let a = Solid::new(first, Operand::First, Some(base))?;
185    let b = Solid::new(second, Operand::Second, Some(base))?;
186    Ok((a, b))
187}
188
189/// The shared volume, clamped to what either solid holds.
190fn shared(a: &Solid, b: &Solid) -> VolumeInterval {
191    let total = pair_sum(&a.faces, &b.faces, a.base);
192    clamp(total, 0.0, a.volume.hi.min(b.volume.hi))
193}
194
195/// `sum over f, g of  s_f s_g * integral over f' and g' of (min(h_f, h_g) - base)`,
196/// unclamped.
197pub(crate) fn pair_sum(first: &[Face], second: &[Face], base: f64) -> Iv {
198    let items: Vec<SpatialItem<u32>> = second
199        .iter()
200        .enumerate()
201        .map(|(i, f)| SpatialItem::new(i as u32, f.shadow()))
202        .collect();
203    let bvh = Bvh::build(items);
204    let mut total = Iv::point(0.0);
205    for f in first {
206        bvh.visit_aabb(&f.shadow(), &mut |j: &u32| {
207            let g = &second[*j as usize];
208            let term = overlap(f, g, base).within(0.0, cap(f, g, base));
209            total = if f.up == g.up {
210                total.add(term)
211            } else {
212                total.sub(term)
213            };
214            ControlFlow::Continue(())
215        });
216    }
217    total
218}
219
220pub(crate) fn clamp(v: Iv, low: f64, high: f64) -> VolumeInterval {
221    VolumeInterval {
222        lower: v.lo.max(low).min(high),
223        upper: v.hi.min(high).max(low),
224    }
225}
226
227pub(crate) fn lowest<M: TriangleMeshView + ?Sized>(mesh: &M) -> f64 {
228    (0..mesh.position_count())
229        .map(|i| mesh.position(i).z)
230        .fold(f64::INFINITY, f64::min)
231}
232
233/// A face that casts a shadow, its corners counter-clockwise in the
234/// plane, and whether its outward side (after orienting the solid) faces
235/// up.
236#[derive(Debug, Clone, Copy)]
237pub(crate) struct Face {
238    p: [Point3; 3],
239    pub(crate) up: bool,
240}
241
242/// No less than any term of `f` and `g`: the smaller shadow's area times
243/// the lower of the two highest corners, above the base.
244fn cap(f: &Face, g: &Face, base: f64) -> f64 {
245    let area = f.area().hi.min(g.area().hi);
246    let top = f.top().min(g.top());
247    Iv::point(area).mul(Iv::point(top).sub(Iv::point(base))).hi
248}
249
250impl Face {
251    fn shadow(&self) -> Aabb {
252        let mut b = Aabb::from_point(Point3::new(self.p[0].x, self.p[0].y, 0.0));
253        for q in &self.p[1..] {
254            b.extend(Point3::new(q.x, q.y, 0.0));
255        }
256        b
257    }
258
259    /// The shadow's area, enclosed.
260    fn area(&self) -> Iv {
261        triangle_area(self.flat(0), self.flat(1), self.flat(2))
262    }
263
264    /// The highest corner.
265    fn top(&self) -> f64 {
266        self.p[0].z.max(self.p[1].z).max(self.p[2].z)
267    }
268
269    fn flat(&self, i: usize) -> Point2 {
270        Point2::new(self.p[i].x, self.p[i].y)
271    }
272
273    /// The edge lines, each from a corner to the next.
274    fn lines(&self) -> [Line; 3] {
275        [0, 1, 2].map(|i| Line {
276            a: self.flat(i),
277            b: self.flat((i + 1) % 3),
278        })
279    }
280
281    /// The plane's normal, enclosed: `(p1 - p0) x (p2 - p0)`.
282    fn normal(&self) -> [Iv; 3] {
283        let d = |i: usize, k: usize| Iv::point(self.p[i][k]).sub(Iv::point(self.p[0][k]));
284        let (u, v) = ([d(1, 0), d(1, 1), d(1, 2)], [d(2, 0), d(2, 1), d(2, 2)]);
285        [
286            u[1].mul(v[2]).sub(u[2].mul(v[1])),
287            u[2].mul(v[0]).sub(u[0].mul(v[2])),
288            u[0].mul(v[1]).sub(u[1].mul(v[0])),
289        ]
290    }
291
292    /// Height of the plane above `(x, y)`, less `base`.
293    fn height(&self, at: [Iv; 2], base: f64) -> Iv {
294        let n = self.normal();
295        let dx = at[0].sub(Iv::point(self.p[0].x));
296        let dy = at[1].sub(Iv::point(self.p[0].y));
297        let lift = n[0].mul(dx).add(n[1].mul(dy)).div(n[2]);
298        Iv::point(self.p[0].z).sub(Iv::point(base)).sub(lift)
299    }
300}
301
302struct Solid {
303    faces: Vec<Face>,
304    volume: Iv,
305    base: f64,
306}
307
308impl Solid {
309    fn new<M: TriangleMeshView + ?Sized>(
310        mesh: &M,
311        operand: Operand,
312        base: Option<f64>,
313    ) -> Result<Self, OverlapError> {
314        if (0..mesh.position_count()).any(|i| !mesh.position(i).is_finite()) {
315            return Err(OverlapError::NonFinite { operand });
316        }
317        if !audit_mesh(mesh, Tolerance::ZERO).is_closed_two_manifold() {
318            return Err(OverlapError::NotClosed { operand });
319        }
320        if let Some(pair) = self_intersections(mesh).first() {
321            return Err(OverlapError::SelfIntersecting {
322                operand,
323                triangles: [pair.first, pair.second],
324            });
325        }
326        let base = base.unwrap_or_else(|| lowest(mesh));
327        let mut faces = faces_as_given(mesh);
328        let mut volume = signed_volume(&faces, base);
329        // Wound inward: turn every face over.
330        if volume.hi < 0.0 {
331            volume = Iv::point(0.0).sub(volume);
332            for f in &mut faces {
333                f.up = !f.up;
334            }
335        } else if volume.lo <= 0.0 {
336            return Err(OverlapError::NoVolume { operand });
337        }
338        Ok(Self {
339            faces,
340            volume,
341            base,
342        })
343    }
344}
345
346/// The faces that cast a shadow, each wound as the mesh gives it.
347pub(crate) fn faces_as_given<M: TriangleMeshView + ?Sized>(mesh: &M) -> Vec<Face> {
348    let mut faces = Vec::new();
349    for t in 0..mesh.triangle_count() {
350        let [a, b, c] = mesh.triangle(t).map(|i| mesh.position(i as usize));
351        let flat = |p: Point3| Point2::new(p.x, p.y);
352        match sign(&Orient {
353            a: flat(a),
354            b: flat(b),
355            c: flat(c),
356        }) {
357            Sign::Positive => faces.push(Face {
358                p: [a, b, c],
359                up: true,
360            }),
361            Sign::Negative => faces.push(Face {
362                p: [a, c, b],
363                up: false,
364            }),
365            _ => {}
366        }
367    }
368    faces
369}
370
371/// The signed volume of the faces: each shadow's area times the mean
372/// height above `base`.
373pub(crate) fn signed_volume(faces: &[Face], base: f64) -> Iv {
374    let mut volume = Iv::point(0.0);
375    for f in faces {
376        let area = f.area();
377        let rise = [0, 1, 2]
378            .map(|i| Iv::point(f.p[i].z).sub(Iv::point(base)))
379            .into_iter()
380            .fold(Iv::point(0.0), Iv::add)
381            .div(Iv::point(3.0));
382        let term = area.mul(rise);
383        volume = if f.up {
384            volume.add(term)
385        } else {
386            volume.sub(term)
387        };
388    }
389    volume
390}
391
392fn triangle_area(a: Point2, b: Point2, c: Point2) -> Iv {
393    let (ux, uy) = (
394        Iv::point(b.x).sub(Iv::point(a.x)),
395        Iv::point(b.y).sub(Iv::point(a.y)),
396    );
397    let (vx, vy) = (
398        Iv::point(c.x).sub(Iv::point(a.x)),
399        Iv::point(c.y).sub(Iv::point(a.y)),
400    );
401    ux.mul(vy).sub(uy.mul(vx)).mul(Iv::point(0.5))
402}
403
404/// The integral of `min(h_f, h_g) - base` over the overlap of the two
405/// shadows.
406fn overlap(f: &Face, g: &Face, base: f64) -> Iv {
407    // The overlap: f's shadow clipped by each of g's edge lines, exactly.
408    let mut poly: Vec<(Vertex, Line)> = (0..3)
409        .map(|i| (Vertex::Input(f.flat(i)), f.lines()[i]))
410        .collect();
411    for line in g.lines() {
412        poly = clip(&poly, line);
413        if poly.len() < 3 {
414            return Iv::point(0.0);
415        }
416    }
417    // Which plane is lower at each vertex, exactly.
418    let signs: Vec<Sign> = poly
419        .iter()
420        .map(|(v, _)| {
421            sign(&Lower {
422                f: *f,
423                g: *g,
424                v: *v,
425            })
426        })
427        .collect();
428    let points: Vec<[Iv; 2]> = poly.iter().map(|(v, _)| v.enclose()).collect();
429    if signs.iter().all(|s| *s == Sign::Zero) {
430        // Coplanar here: one plane over the whole overlap.
431        return integral(&points, f, base);
432    }
433    let rise = |p: [Iv; 2]| f.height(p, base).sub(g.height(p, base));
434    let below = split(&points, &signs, Sign::Negative, &rise);
435    let above = split(&points, &signs, Sign::Positive, &rise);
436    integral(&below, f, base).add(integral(&above, g, base))
437}
438
439/// The part of a convex polygon where the sign is `keep` or zero, cut
440/// where the plane difference changes sign.
441fn split(
442    points: &[[Iv; 2]],
443    signs: &[Sign],
444    keep: Sign,
445    rise: &dyn Fn([Iv; 2]) -> Iv,
446) -> Vec<[Iv; 2]> {
447    let n = points.len();
448    let mut out = Vec::new();
449    for i in 0..n {
450        let j = (i + 1) % n;
451        let (si, sj) = (signs[i], signs[j]);
452        if si == keep || si == Sign::Zero {
453            out.push(points[i]);
454        }
455        let opposite = matches!(
456            (si, sj),
457            (Sign::Negative, Sign::Positive) | (Sign::Positive, Sign::Negative)
458        );
459        if opposite {
460            // Where along the edge the difference is zero: in (0, 1).
461            let (ri, rj) = (rise(points[i]), rise(points[j]));
462            let t = ri.div(ri.sub(rj)).within(0.0, 1.0);
463            let at = |k: usize| points[i][k].add(t.mul(points[j][k].sub(points[i][k])));
464            out.push([at(0), at(1)]);
465        }
466    }
467    out
468}
469
470/// The integral of the face's height above the base over a convex
471/// polygon: a fan of triangles, each its area times its mean height.
472fn integral(points: &[[Iv; 2]], face: &Face, base: f64) -> Iv {
473    if points.len() < 3 {
474        return Iv::point(0.0);
475    }
476    let heights: Vec<Iv> = points.iter().map(|p| face.height(*p, base)).collect();
477    let mut total = Iv::point(0.0);
478    for i in 1..points.len() - 1 {
479        let (a, b, c) = (points[0], points[i], points[i + 1]);
480        let area = b[0]
481            .sub(a[0])
482            .mul(c[1].sub(a[1]))
483            .sub(b[1].sub(a[1]).mul(c[0].sub(a[0])))
484            .mul(Iv::point(0.5));
485        let mean = heights[0]
486            .add(heights[i])
487            .add(heights[i + 1])
488            .div(Iv::point(3.0));
489        total = total.add(area.mul(mean));
490    }
491    total
492}
493
494/// A line through two input points, directed from `a` to `b`.
495#[derive(Debug, Clone, Copy, PartialEq)]
496struct Line {
497    a: Point2,
498    b: Point2,
499}
500
501/// A vertex of an overlap: an input corner, or where two input lines
502/// cross.
503#[derive(Debug, Clone, Copy)]
504enum Vertex {
505    Input(Point2),
506    Cross(Line, Line),
507}
508
509impl Vertex {
510    /// Homogeneous coordinates `(X, Y, W)` in the arithmetic `T`: the
511    /// point is `(X / W, Y / W)`, with `W` the lines' cross product for a
512    /// crossing.
513    fn homogeneous<T: Arith>(&self) -> [T; 3] {
514        let f = T::from_f64;
515        match *self {
516            Vertex::Input(p) => [f(p.x), f(p.y), f(1.0)],
517            Vertex::Cross(l1, l2) => {
518                let d1 = [f(l1.b.x).sub(&f(l1.a.x)), f(l1.b.y).sub(&f(l1.a.y))];
519                let d2 = [f(l2.b.x).sub(&f(l2.a.x)), f(l2.b.y).sub(&f(l2.a.y))];
520                let w = d1[0].mul(&d2[1]).sub(&d1[1].mul(&d2[0]));
521                let e = [f(l2.a.x).sub(&f(l1.a.x)), f(l2.a.y).sub(&f(l1.a.y))];
522                let n = e[0].mul(&d2[1]).sub(&e[1].mul(&d2[0]));
523                [
524                    f(l1.a.x).mul(&w).add(&d1[0].mul(&n)),
525                    f(l1.a.y).mul(&w).add(&d1[1].mul(&n)),
526                    w,
527                ]
528            }
529        }
530    }
531
532    fn enclose(&self) -> [Iv; 2] {
533        match *self {
534            Vertex::Input(p) => [Iv::point(p.x), Iv::point(p.y)],
535            Vertex::Cross(l1, l2) => {
536                let d1 = [
537                    Iv::point(l1.b.x).sub(Iv::point(l1.a.x)),
538                    Iv::point(l1.b.y).sub(Iv::point(l1.a.y)),
539                ];
540                let d2 = [
541                    Iv::point(l2.b.x).sub(Iv::point(l2.a.x)),
542                    Iv::point(l2.b.y).sub(Iv::point(l2.a.y)),
543                ];
544                let w = d1[0].mul(d2[1]).sub(d1[1].mul(d2[0]));
545                let e = [
546                    Iv::point(l2.a.x).sub(Iv::point(l1.a.x)),
547                    Iv::point(l2.a.y).sub(Iv::point(l1.a.y)),
548                ];
549                let t = e[0].mul(d2[1]).sub(e[1].mul(d2[0])).div(w);
550                [
551                    Iv::point(l1.a.x).add(d1[0].mul(t)),
552                    Iv::point(l1.a.y).add(d1[1].mul(t)),
553                ]
554            }
555        }
556    }
557}
558
559/// Sutherland-Hodgman against the closed left side of `line`, keeping
560/// the line under each edge so every new vertex is a crossing of two
561/// input lines.
562fn clip(poly: &[(Vertex, Line)], line: Line) -> Vec<(Vertex, Line)> {
563    let n = poly.len();
564    let signs: Vec<Sign> = poly
565        .iter()
566        .map(|(v, _)| sign(&SideOf { line, v: *v }))
567        .collect();
568    let mut out = Vec::with_capacity(n + 1);
569    for i in 0..n {
570        let (v, edge) = poly[i];
571        let (si, sj) = (signs[i], signs[(i + 1) % n]);
572        if si != Sign::Negative {
573            // From a vertex on the line to the next one out, the kept
574            // boundary runs along the line to where it comes back in.
575            let along = si == Sign::Zero && sj == Sign::Negative;
576            out.push((v, if along { line } else { edge }));
577        }
578        match (si, sj) {
579            (Sign::Positive, Sign::Negative) => out.push((Vertex::Cross(edge, line), line)),
580            (Sign::Negative, Sign::Positive) => out.push((Vertex::Cross(edge, line), edge)),
581            _ => {}
582        }
583    }
584    out
585}
586
587/// The inputs are finite, so the exact tier always decides.
588fn sign<E: SignExpr>(e: &E) -> Sign {
589    certify(e).unwrap_or(Sign::Zero)
590}
591
592struct Orient {
593    a: Point2,
594    b: Point2,
595    c: Point2,
596}
597
598impl SignExpr for Orient {
599    fn sign_in<T: Arith>(&self) -> Option<Sign> {
600        SideOf {
601            line: Line {
602                a: self.a,
603                b: self.b,
604            },
605            v: Vertex::Input(self.c),
606        }
607        .sign_in::<T>()
608    }
609}
610
611/// Which side of `line` a vertex lies on: positive to the left.
612struct SideOf {
613    line: Line,
614    v: Vertex,
615}
616
617impl SignExpr for SideOf {
618    fn sign_in<T: Arith>(&self) -> Option<Sign> {
619        let f = T::from_f64;
620        let [x, y, w] = self.v.homogeneous::<T>();
621        let (a, b) = (self.line.a, self.line.b);
622        let (dx, dy) = (f(b.x).sub(&f(a.x)), f(b.y).sub(&f(a.y)));
623        let rx = x.sub(&f(a.x).mul(&w));
624        let ry = y.sub(&f(a.y).mul(&w));
625        let s = dx.mul(&ry).sub(&dy.mul(&rx)).sign()?;
626        Some(times(s, w.sign()?))
627    }
628}
629
630/// The sign of `h_f - h_g` at a vertex: negative where `f`'s plane is the
631/// lower. Both shadows are counter-clockwise, so both planes' normals
632/// point up and multiplying through by them keeps the sign.
633struct Lower {
634    f: Face,
635    g: Face,
636    v: Vertex,
637}
638
639impl SignExpr for Lower {
640    fn sign_in<T: Arith>(&self) -> Option<Sign> {
641        let f = T::from_f64;
642        let [x, y, w] = self.v.homogeneous::<T>();
643        let normal = |face: &Face| {
644            let d = |i: usize, k: usize| f(face.p[i][k]).sub(&f(face.p[0][k]));
645            let (u, v) = ([d(1, 0), d(1, 1), d(1, 2)], [d(2, 0), d(2, 1), d(2, 2)]);
646            [
647                u[1].mul(&v[2]).sub(&u[2].mul(&v[1])),
648                u[2].mul(&v[0]).sub(&u[0].mul(&v[2])),
649                u[0].mul(&v[1]).sub(&u[1].mul(&v[0])),
650            ]
651        };
652        let (nf, ng) = (normal(&self.f), normal(&self.g));
653        // W n_z h(x, y) = W n_z p0.z - n_x (X - W p0.x) - n_y (Y - W p0.y).
654        let lift = |face: &Face, n: &[T; 3]| {
655            let p = face.p[0];
656            w.mul(&n[2])
657                .mul(&f(p.z))
658                .sub(&n[0].mul(&x.sub(&w.mul(&f(p.x)))))
659                .sub(&n[1].mul(&y.sub(&w.mul(&f(p.y)))))
660        };
661        let difference = lift(&self.f, &nf)
662            .mul(&ng[2])
663            .sub(&lift(&self.g, &ng).mul(&nf[2]));
664        Some(times(difference.sign()?, w.sign()?))
665    }
666}
667
668fn times(a: Sign, b: Sign) -> Sign {
669    match (a, b) {
670        (Sign::Zero, _) | (_, Sign::Zero) => Sign::Zero,
671        (x, y) if x == y => Sign::Positive,
672        _ => Sign::Negative,
673    }
674}
675
676/// An outward-rounded interval: every operation rounds to nearest and
677/// steps each bound one float outward, which covers the rounding error.
678#[derive(Debug, Clone, Copy, PartialEq)]
679pub(crate) struct Iv {
680    pub(crate) lo: f64,
681    pub(crate) hi: f64,
682}
683
684impl Iv {
685    pub(crate) fn point(v: f64) -> Self {
686        Self { lo: v, hi: v }
687    }
688
689    fn new(lo: f64, hi: f64) -> Self {
690        Self { lo, hi }
691    }
692
693    fn outward(lo: f64, hi: f64) -> Self {
694        if lo.is_nan() || hi.is_nan() {
695            return Self::new(f64::NEG_INFINITY, f64::INFINITY);
696        }
697        Self::new(lo.next_down(), hi.next_up())
698    }
699
700    pub(crate) fn add(self, o: Self) -> Self {
701        Self::outward(self.lo + o.lo, self.hi + o.hi)
702    }
703
704    pub(crate) fn sub(self, o: Self) -> Self {
705        Self::outward(self.lo - o.hi, self.hi - o.lo)
706    }
707
708    pub(crate) fn mul(self, o: Self) -> Self {
709        let p = [
710            self.lo * o.lo,
711            self.lo * o.hi,
712            self.hi * o.lo,
713            self.hi * o.hi,
714        ];
715        let lo = p.iter().copied().fold(f64::INFINITY, f64::min);
716        let hi = p.iter().copied().fold(f64::NEG_INFINITY, f64::max);
717        Self::outward(lo, hi)
718    }
719
720    /// Quotient; the whole line if the divisor may be zero.
721    fn div(self, o: Self) -> Self {
722        if o.lo <= 0.0 && o.hi >= 0.0 {
723            return Self::new(f64::NEG_INFINITY, f64::INFINITY);
724        }
725        let q = [
726            self.lo / o.lo,
727            self.lo / o.hi,
728            self.hi / o.lo,
729            self.hi / o.hi,
730        ];
731        let lo = q.iter().copied().fold(f64::INFINITY, f64::min);
732        let hi = q.iter().copied().fold(f64::NEG_INFINITY, f64::max);
733        Self::outward(lo, hi)
734    }
735
736    /// Intersected with `[low, high]`, where the value is known to lie.
737    fn within(self, low: f64, high: f64) -> Self {
738        Self::new(self.lo.max(low), self.hi.min(high))
739    }
740}
axiolid_inspect/overlap.rs

axiolid_inspect/
overlap.rs
