ogeom_algo/
certify.rs

1//! Proving that a face lies within a region bounded by planes.
2//!
3//! A [`HalfSpaces`] is a tree of half-spaces under maxima and minima: a convex
4//! corner is the maximum of its faces' heights, a concave edge the minimum
5//! of its two, and so on. [`certify_face_within`] proves that every point
6//! of a face stands no farther past the region than a stated fit, or finds
7//! a point that does.
8//!
9//! The proof rests on the convex hull. Every surface kind converts exactly
10//! to a B-spline over any rectangle of its own parameters
11//! ([`SurfaceGeometry::to_bspline`] of a [`TrimmedSurface`]), and a B-spline
12//! whose weights are all positive lies in the convex hull of its control
13//! points. A plane's height is affine, so its greatest over the hull is its
14//! greatest at a control point. The maximum of several pieces is bounded by
15//! the maximum of their bounds, and the minimum by the minimum of theirs. A
16//! rectangle whose bound is within the fit is cleared; one wholly outside
17//! the face's trim is passed over; the rest is split in four and tried
18//! again, down to [`DEPTH`] halvings each way and [`BUDGET`] rectangles in
19//! all.
20//!
21//! The trim is read from the face's rings in its chart, each drawn within a
22//! chord of its edges in space through the surface's tangent scales
23//! ([`ogeom_mesh::face_boundary_in_chart_met`]). A rectangle is passed over
24//! only where its centre lies outside the rings and every ring stands
25//! farther from that centre than the rectangle's half diagonal and a slack
26//! of [`SLACK`] chords read in the chart at that centre. That reading of
27//! the trim is the one input the proof takes from the mesher.
28
29use ogeom_core::{OgeomResult, Tolerances};
30use ogeom_geom::{BSplineSurface, Surface as _, SurfaceGeometry, TrimmedSurface};
31use ogeom_math::{Point, Point2, Vector};
32use ogeom_topo::{Model, Shape};
33
34/// A region of space bounded by planes, as a tree of half-spaces.
35#[derive(Debug, Clone, PartialEq)]
36pub enum HalfSpaces {
37    /// The points on the side of the plane through `at` that `normal`
38    /// points away from: a point's height is `normal · (p - at)`, with
39    /// `normal` of unit length.
40    Below {
41        /// The plane's outward unit normal.
42        normal: Vector,
43        /// A point on the plane.
44        at: Point,
45    },
46    /// The points inside every piece: a height the greatest of the
47    /// pieces'.
48    Every(Vec<HalfSpaces>),
49    /// The points inside any piece: a height the least of the pieces'.
50    Any(Vec<HalfSpaces>),
51}
52
53impl HalfSpaces {
54    /// How far `p` stands past the region (negative inside), and the
55    /// outward normal of the plane that decides it.
56    #[must_use]
57    pub fn height(&self, p: Point) -> (f64, Vector) {
58        match self {
59            Self::Below { normal, at } => (normal.dot(p - *at), *normal),
60            Self::Every(pieces) => pieces.iter().map(|r| r.height(p)).fold(
61                (f64::NEG_INFINITY, Vector::ZERO),
62                |a, b| {
63                    if b.0 > a.0 { b } else { a }
64                },
65            ),
66            Self::Any(pieces) => {
67                pieces
68                    .iter()
69                    .map(|r| r.height(p))
70                    .fold(
71                        (f64::INFINITY, Vector::ZERO),
72                        |a, b| {
73                            if b.0 < a.0 { b } else { a }
74                        },
75                    )
76            }
77        }
78    }
79
80    /// A bound on the height over the convex hull of `hull`: exact for a
81    /// plane and for the greatest of pieces, at least the least of pieces'
82    /// heights for the least.
83    #[must_use]
84    pub fn bound_over(&self, hull: &[Point]) -> f64 {
85        match self {
86            Self::Below { normal, at } => hull
87                .iter()
88                .map(|p| normal.dot(*p - *at))
89                .fold(f64::NEG_INFINITY, f64::max),
90            Self::Every(pieces) => pieces
91                .iter()
92                .map(|r| r.bound_over(hull))
93                .fold(f64::NEG_INFINITY, f64::max),
94            Self::Any(pieces) => pieces
95                .iter()
96                .map(|r| r.bound_over(hull))
97                .fold(f64::INFINITY, f64::min),
98        }
99    }
100}
101
102/// What kept a face from being proved within its region.
103#[derive(Debug, Clone, Copy, PartialEq)]
104pub struct Uncertified {
105    /// How far past the region: the height at `at` where `found`, else the
106    /// bound over the rectangle the proof could not clear.
107    pub past: f64,
108    /// A point of the face standing `past` outside, or the centre of the
109    /// rectangle the proof could not clear.
110    pub at: Point,
111    /// Whether `at` is a point of the face found standing past the fit, as
112    /// against one the proof ran out of halvings over.
113    pub found: bool,
114}
115
116/// How many times the face's parameter rectangle is halved each way before
117/// a part of it no bound clears is taken as standing outside.
118pub const DEPTH: u32 = 20;
119
120/// How many rectangles a proof may try before what it has not cleared is
121/// taken as standing outside.
122pub const BUDGET: usize = 1 << 20;
123
124/// The chords of the trim's drawing allowed between a rectangle and the
125/// rings before it is passed over as outside the trim.
126pub const SLACK: f64 = 2.0;
127
128/// Whether every point of `face` stands within `fit` of `region`: `Ok(())`
129/// when proved, else what kept it from being proved (see the module
130/// documentation for the argument and its one input from the mesher).
131///
132/// # Errors
133///
134/// As reading the face's surface and its rings, or converting a rectangle
135/// of the surface to a B-spline.
136pub fn certify_face_within(
137    model: &Model,
138    face: &Shape,
139    region: &HalfSpaces,
140    fit: f64,
141    tol: Tolerances,
142) -> OgeomResult<Result<(), Uncertified>> {
143    let Some(data) = model.node(face).and_then(|n| n.data().as_face()) else {
144        ogeom_core::ogeom_bail!(Dangling, "face is not in this model");
145    };
146    let Some(surface) = model.geometry().surface(data.surface).cloned() else {
147        ogeom_core::ogeom_bail!(Dangling, "face refers to a surface not in this model");
148    };
149    let place = face.transform(model.datums())?;
150    let chord = fit / 4.0;
151    let deflection = ogeom_mesh::Deflection {
152        chord,
153        ..ogeom_mesh::Deflection::default()
154    };
155    let Some(rings) = ogeom_mesh::face_boundary_in_chart_met(model, face, deflection, tol)? else {
156        // A trim drawn short of its chord says nothing the proof can use.
157        let ((u, _), (v, _)) = surface.domain();
158        return Ok(Err(Uncertified {
159            past: f64::INFINITY,
160            at: place.apply(surface.point_at(u, v, tol)?),
161            found: false,
162        }));
163    };
164    let Some(trim) = Trim::of(&rings) else {
165        return Ok(Ok(()));
166    };
167    // A B-spline's pieces are cut from their parent's, keeping its
168    // parameters; any other surface's converted afresh over each rectangle.
169    let spline = match &surface {
170        SurfaceGeometry::BSpline(s) => {
171            let ((u0, u1), (v0, v1)) = s.domain();
172            let (lo, hi) = (trim.lo, trim.hi);
173            Some(s.segment(
174                (lo.x.max(u0), hi.x.min(u1)),
175                (lo.y.max(v0), hi.y.min(v1)),
176                tol,
177            )?)
178        }
179        _ => None,
180    };
181    let start = match &spline {
182        Some(s) => {
183            let ((u0, u1), (v0, v1)) = s.domain();
184            (Point2::new(u0, v0), Point2::new(u1, v1))
185        }
186        None => (trim.lo, trim.hi),
187    };
188    let mut pending: Vec<(Point2, Point2, u32, Option<BSplineSurface>)> =
189        vec![(start.0, start.1, 0, spline)];
190    let mut tried = 0usize;
191    while let Some((a, b, depth, parent)) = pending.pop() {
192        tried += 1;
193        let centre = Point2::new(f64::midpoint(a.x, b.x), f64::midpoint(a.y, b.y));
194        let inside = trim.inside(centre);
195        if !inside {
196            // The chart's slack at the centre: chords of space read back
197            // through the tangent scales there.
198            let (du, dv) = surface.d1_at(centre.x, centre.y, tol)?;
199            let scale = du.magnitude().min(dv.magnitude());
200            let slack = if scale > 0.0 {
201                SLACK * chord / scale
202            } else {
203                f64::INFINITY
204            };
205            let half = (b.x - a.x).hypot(b.y - a.y) * 0.5;
206            if !trim.near(centre, half + slack) {
207                continue;
208            }
209        }
210        let piece = match &parent {
211            Some(p) => p.segment((a.x, b.x), (a.y, b.y), tol)?,
212            None => SurfaceGeometry::Trimmed(Box::new(TrimmedSurface::new(
213                surface.clone(),
214                (a.x, b.x),
215                (a.y, b.y),
216                tol,
217            )?))
218            .to_bspline(tol)?,
219        };
220        let net = piece.grid().points();
221        let bound = if net.iter().all(|w| w.weight > 0.0) {
222            let hull: Vec<Point> = net.iter().map(|w| place.apply(w.point())).collect();
223            region.bound_over(&hull)
224        } else {
225            f64::INFINITY
226        };
227        if bound <= fit {
228            continue;
229        }
230        let at = place.apply(surface.point_at(centre.x, centre.y, tol)?);
231        if inside {
232            let (past, _) = region.height(at);
233            if past > fit {
234                return Ok(Err(Uncertified {
235                    past,
236                    at,
237                    found: true,
238                }));
239            }
240        }
241        if depth >= DEPTH || tried >= BUDGET {
242            return Ok(Err(Uncertified {
243                past: bound,
244                at,
245                found: false,
246            }));
247        }
248        let carried = parent.is_some().then_some(piece);
249        for (x0, x1) in [(a.x, centre.x), (centre.x, b.x)] {
250            for (y0, y1) in [(a.y, centre.y), (centre.y, b.y)] {
251                pending.push((
252                    Point2::new(x0, y0),
253                    Point2::new(x1, y1),
254                    depth + 1,
255                    carried.clone(),
256                ));
257            }
258        }
259    }
260    Ok(Ok(()))
261}
262
263/// A face's trim rings in its chart, bucketed: their segments on a grid of
264/// cells for the nearest within a reach, and by bands of height for the
265/// crossings a ray to `+x` makes.
266struct Trim {
267    segments: Vec<(Point2, Point2)>,
268    lo: Point2,
269    hi: Point2,
270    step: (f64, f64),
271    cells: Vec<Vec<usize>>,
272    bands: Vec<Vec<usize>>,
273}
274
275/// The cells of a [`Trim`] each way.
276const CELLS: usize = 64;
277
278impl Trim {
279    /// The rings bucketed; `None` where they hold no segment.
280    fn of(rings: &[Vec<Point2>]) -> Option<Self> {
281        let segments: Vec<(Point2, Point2)> = rings
282            .iter()
283            .flat_map(|ring| (0..ring.len()).map(move |k| (ring[k], ring[(k + 1) % ring.len()])))
284            .collect();
285        if segments.is_empty() {
286            return None;
287        }
288        let (mut lo, mut hi) = (
289            Point2::new(f64::INFINITY, f64::INFINITY),
290            Point2::new(f64::NEG_INFINITY, f64::NEG_INFINITY),
291        );
292        for (a, _) in &segments {
293            lo = Point2::new(lo.x.min(a.x), lo.y.min(a.y));
294            hi = Point2::new(hi.x.max(a.x), hi.y.max(a.y));
295        }
296        #[expect(clippy::cast_precision_loss, reason = "a cell count of sixty-four")]
297        let step = (
298            ((hi.x - lo.x) / CELLS as f64).max(f64::MIN_POSITIVE),
299            ((hi.y - lo.y) / CELLS as f64).max(f64::MIN_POSITIVE),
300        );
301        let mut trim = Self {
302            segments,
303            lo,
304            hi,
305            step,
306            cells: vec![Vec::new(); CELLS * CELLS],
307            bands: vec![Vec::new(); CELLS],
308        };
309        for k in 0..trim.segments.len() {
310            let (a, b) = trim.segments[k];
311            let (i0, j0) = trim.cell(Point2::new(a.x.min(b.x), a.y.min(b.y)));
312            let (i1, j1) = trim.cell(Point2::new(a.x.max(b.x), a.y.max(b.y)));
313            for i in i0..=i1 {
314                for j in j0..=j1 {
315                    trim.cells[i * CELLS + j].push(k);
316                }
317            }
318            for band in &mut trim.bands[j0..=j1] {
319                band.push(k);
320            }
321        }
322        Some(trim)
323    }
324
325    /// The cell holding `p`, clamped to the grid.
326    fn cell(&self, p: Point2) -> (usize, usize) {
327        #[expect(
328            clippy::cast_possible_truncation,
329            clippy::cast_sign_loss,
330            reason = "clamped to the grid's sixty-four cells"
331        )]
332        let at = |x: f64, lo: f64, step: f64| {
333            (((x - lo) / step).floor().max(0.0) as usize).min(CELLS - 1)
334        };
335        (
336            at(p.x, self.lo.x, self.step.0),
337            at(p.y, self.lo.y, self.step.1),
338        )
339    }
340
341    /// Whether a segment of the rings comes within `reach` of `p`.
342    fn near(&self, p: Point2, reach: f64) -> bool {
343        if !reach.is_finite() {
344            return true;
345        }
346        if p.x + reach < self.lo.x
347            || p.x - reach > self.hi.x
348            || p.y + reach < self.lo.y
349            || p.y - reach > self.hi.y
350        {
351            return false;
352        }
353        let (i0, j0) = self.cell(Point2::new(p.x - reach, p.y - reach));
354        let (i1, j1) = self.cell(Point2::new(p.x + reach, p.y + reach));
355        (i0..=i1).any(|i| {
356            (j0..=j1).any(|j| {
357                self.cells[i * CELLS + j]
358                    .iter()
359                    .any(|&k| distance_to(p, self.segments[k]) <= reach)
360            })
361        })
362    }
363
364    /// Whether `p` lies inside the rings, by the parity of the crossings
365    /// its ray to `+x` makes: only segments of its band can cross it.
366    fn inside(&self, p: Point2) -> bool {
367        if p.y < self.lo.y || p.y > self.hi.y {
368            return false;
369        }
370        let (_, j) = self.cell(p);
371        let mut crossings = 0usize;
372        for &k in &self.bands[j] {
373            let (a, b) = self.segments[k];
374            if (a.y > p.y) != (b.y > p.y) {
375                let x = a.x + (p.y - a.y) / (b.y - a.y) * (b.x - a.x);
376                if x > p.x {
377                    crossings += 1;
378                }
379            }
380        }
381        crossings % 2 == 1
382    }
383}
384
385/// The distance from `p` to the segment `(a, b)`.
386fn distance_to(p: Point2, (a, b): (Point2, Point2)) -> f64 {
387    let (ux, uy) = (b.x - a.x, b.y - a.y);
388    let length = ux.mul_add(ux, uy * uy);
389    let t = if length > 0.0 {
390        (((p.x - a.x) * ux + (p.y - a.y) * uy) / length).clamp(0.0, 1.0)
391    } else {
392        0.0
393    };
394    (p.x - (a.x + ux * t)).hypot(p.y - (a.y + uy * t))
395}
396
397#[cfg(test)]
398mod tests {
399    #![allow(clippy::unwrap_used, reason = "test code")]
400    use super::*;
401    use ogeom_geom::BSplineSurface;
402    use ogeom_math::{ControlGrid, Frame, KnotVector};
403    use ogeom_topo::{ShapeType, explore_unique};
404
405    const T: Tolerances = Tolerances::millimetres();
406    const FIT: f64 = 2e-4;
407
408    fn below(normal: Vector, at: Point) -> HalfSpaces {
409        HalfSpaces::Below { normal, at }
410    }
411
412    /// The square prism `|x|, |y| <= half`, as the maximum of its four sides.
413    fn square(half: f64) -> HalfSpaces {
414        let side =
415            |x: f64, y: f64| below(Vector::new(x, y, 0.0), Point::new(x * half, y * half, 0.0));
416        HalfSpaces::Every(vec![
417            side(1.0, 0.0),
418            side(-1.0, 0.0),
419            side(0.0, 1.0),
420            side(0.0, -1.0),
421        ])
422    }
423
424    /// The faces of `solid` on a surface of `kind`.
425    fn faces_of(model: &Model, solid: &Shape, kind: ogeom_geom::SurfaceKind) -> Vec<Shape> {
426        explore_unique(model, solid, ShapeType::Face)
427            .unwrap()
428            .into_iter()
429            .filter(|f| {
430                let data = model.node(f).unwrap().data().as_face().unwrap();
431                model.geometry().surface(data.surface).unwrap().kind() == kind
432            })
433            .collect()
434    }
435
436    /// A unit cylinder's side touches the square prism round it along four
437    /// lines and lies inside it: proved, the tangency cleared by the arcs'
438    /// middle control points standing on the sides. Drawn in by a tenth,
439    /// the sides cut it, and a point standing out is found.
440    #[test]
441    fn a_cylinder_inside_its_square_is_proved_and_one_poking_out_is_found() {
442        let mut m = Model::new();
443        let solid = crate::make_cylinder(&mut m, Frame::WORLD, 1.0, 3.0, T)
444            .unwrap()
445            .shape;
446        let side = faces_of(&m, &solid, ogeom_geom::SurfaceKind::Cylinder);
447        assert_eq!(side.len(), 1);
448        assert_eq!(
449            certify_face_within(&m, &side[0], &square(1.0), FIT, T).unwrap(),
450            Ok(())
451        );
452        let out = certify_face_within(&m, &side[0], &square(0.9), FIT, T)
453            .unwrap()
454            .unwrap_err();
455        assert!(out.found, "{out:?}");
456        assert!(out.past > FIT && out.past <= 0.1 + 1e-12, "{out:?}");
457    }
458
459    /// A sphere of radius 2 inside the cube it touches, the cube as the
460    /// maximum of six planes: proved. With one side drawn in by 0.05, the
461    /// point found stands past it, by no more than 0.05.
462    #[test]
463    fn a_sphere_in_its_cube_is_proved_and_past_a_drawn_in_side_is_found() {
464        let mut m = Model::new();
465        let solid = crate::make_sphere(&mut m, Frame::WORLD, 2.0, T)
466            .unwrap()
467            .shape;
468        let faces = faces_of(&m, &solid, ogeom_geom::SurfaceKind::Sphere);
469        let side = |n: Vector, d: f64| below(n, Point::ORIGIN + n * d);
470        let axes = [
471            Vector::new(1.0, 0.0, 0.0),
472            Vector::new(-1.0, 0.0, 0.0),
473            Vector::new(0.0, 1.0, 0.0),
474            Vector::new(0.0, -1.0, 0.0),
475            Vector::new(0.0, 0.0, 1.0),
476            Vector::new(0.0, 0.0, -1.0),
477        ];
478        let cube = HalfSpaces::Every(axes.iter().map(|n| side(*n, 2.0)).collect());
479        for f in &faces {
480            assert_eq!(certify_face_within(&m, f, &cube, FIT, T).unwrap(), Ok(()));
481        }
482        let mut pieces: Vec<HalfSpaces> = axes.iter().map(|n| side(*n, 2.0)).collect();
483        pieces[4] = side(axes[4], 1.95);
484        let drawn = HalfSpaces::Every(pieces);
485        let worst = faces
486            .iter()
487            .filter_map(|f| certify_face_within(&m, f, &drawn, FIT, T).unwrap().err())
488            .collect::<Vec<_>>();
489        assert!(!worst.is_empty());
490        assert!(
491            worst
492                .iter()
493                .all(|u| u.found && u.past > FIT && u.past <= 0.05 + 1e-12)
494        );
495    }
496
497    /// The soundness hole a guard reading only vertices leaves: a
498    /// biquadratic patch over the unit square whose corners and edges lie
499    /// on the plane `z = 0` but whose middle control point stands 1 above
500    /// it. Every vertex is on the plane; the patch stands a quarter past it
501    /// at its middle (the middle basis function is a half squared there),
502    /// and that is what is found, to the bit.
503    #[test]
504    fn a_patch_bulging_past_its_plane_between_its_vertices_is_found() {
505        let mut m = Model::new();
506        let mut points = Vec::new();
507        for i in 0..3 {
508            for j in 0..3 {
509                let z = if i == 1 && j == 1 { 1.0 } else { 0.0 };
510                points.push(Point::new(f64::from(i) * 0.5, f64::from(j) * 0.5, z));
511            }
512        }
513        let knots = || KnotVector::new(vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0], 2).unwrap();
514        let grid = ControlGrid::new(points, 3, 3).unwrap();
515        let surface = BSplineSurface::new(knots(), knots(), &grid, T).unwrap();
516        let face = crate::make_natural_face(&mut m, surface.into())
517            .unwrap()
518            .shape;
519        let plane = below(Vector::new(0.0, 0.0, 1.0), Point::ORIGIN);
520        let out = certify_face_within(&m, &face, &plane, FIT, T)
521            .unwrap()
522            .unwrap_err();
523        assert!(out.found, "{out:?}");
524        assert_eq!(out.past, 0.25, "{out:?}");
525        // Raised by the quarter, the plane holds it.
526        let raised = below(Vector::new(0.0, 0.0, 1.0), Point::new(0.0, 0.0, 0.25));
527        assert_eq!(
528            certify_face_within(&m, &face, &raised, FIT, T).unwrap(),
529            Ok(())
530        );
531    }
532
533    /// The least of pieces: a point is inside where it is inside any. The
534    /// unit cylinder's side lies inside the union of the half-spaces
535    /// `x <= 1` and `x >= -1` everywhere, so proved; and inside the union of
536    /// `x <= -0.5` and `x >= 0.5` nowhere near `x = 0`, so a point is found.
537    #[test]
538    fn the_least_of_pieces_proves_a_union_and_finds_what_falls_between() {
539        let mut m = Model::new();
540        let solid = crate::make_cylinder(&mut m, Frame::WORLD, 1.0, 3.0, T)
541            .unwrap()
542            .shape;
543        let side = &faces_of(&m, &solid, ogeom_geom::SurfaceKind::Cylinder)[0];
544        let x = Vector::new(1.0, 0.0, 0.0);
545        let both = HalfSpaces::Any(vec![
546            below(x, Point::new(1.0, 0.0, 0.0)),
547            below(-x, Point::new(-1.0, 0.0, 0.0)),
548        ]);
549        assert_eq!(
550            certify_face_within(&m, side, &both, FIT, T).unwrap(),
551            Ok(())
552        );
553        let apart = HalfSpaces::Any(vec![
554            below(x, Point::new(-0.5, 0.0, 0.0)),
555            below(-x, Point::new(0.5, 0.0, 0.0)),
556        ]);
557        let out = certify_face_within(&m, side, &apart, FIT, T)
558            .unwrap()
559            .unwrap_err();
560        assert!(
561            out.found && out.past > FIT && out.past <= 0.5 + 1e-12,
562            "{out:?}"
563        );
564    }
565}
ogeom_algo/certify.rs

ogeom_algo/
certify.rs
