ogeom_algo/
tight.rs

1//! The smallest axis-aligned box that holds a shape: how big a body is.
2//!
3//! [`shape_bounds`](crate::shape_bounds) promises to hold everything and
4//! so keeps each surface's carrier bound, which for a surface of
5//! revolution or a trimmed cone can stand well clear of the solid. Here
6//! each extreme is found where it is: on an edge, by a one-dimensional
7//! search along its curve, or inside a face, at a point of the exact
8//! surface found from the face's own mesh and held to the face's chart.
9
10use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
11use ogeom_geom::{Curve3d as _, Surface as _};
12use ogeom_math::{Aabb, Point, Point2, Vector, solve};
13use ogeom_mesh::Deflection;
14use ogeom_topo::{EdgeRepr, Model, NodeData, Shape, ShapeType, explore_unique};
15
16const AXES: [Vector; 3] = [Vector::X, Vector::Y, Vector::Z];
17
18/// The smallest axis-aligned box holding `shape`, to within `tol`: each of
19/// its six sides where the shape actually reaches, on an edge, at a vertex
20/// or inside a face.
21///
22/// # Errors
23///
24/// [`OgeomError::Dangling`](ogeom_core::OgeomError::Dangling) if a handle
25/// does not resolve; [`OgeomError::Construction`](ogeom_core::OgeomError::Construction)
26/// if the shape holds nothing to bound, or as the face meshes report.
27pub fn tight_bounds(model: &Model, shape: &Shape, tol: Tolerances) -> OgeomResult<Aabb> {
28    let mut points: Vec<Point> = Vec::new();
29
30    for vertex in explore_unique(model, shape, ShapeType::Vertex)? {
31        let Some(data) = model.node(&vertex).and_then(|n| n.data().as_vertex()) else {
32            continue;
33        };
34        points.push(vertex.transform(model.datums())?.apply(data.point));
35    }
36
37    // Along every edge, each coordinate's least and greatest, the search
38    // bracketed by the samples that led.
39    for edge in explore_unique(model, shape, ShapeType::Edge)? {
40        let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
41            continue;
42        };
43        let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
44            continue;
45        };
46        let Some(geometry) = model.geometry().curve(*curve) else {
47            ogeom_bail!(Dangling, "curve is not in this model");
48        };
49        let placement = edge.transform(model.datums())?;
50        let at = |t: f64| -> OgeomResult<Point> { Ok(placement.apply(geometry.point_at(t, tol)?)) };
51        const N: usize = 64;
52        #[allow(clippy::cast_precision_loss)]
53        let step = (range.1 - range.0) / N as f64;
54        let samples: Vec<(f64, Point)> = (0..=N)
55            .map(|i| {
56                #[allow(clippy::cast_precision_loss)]
57                let t = range.0 + step * i as f64;
58                at(t).map(|p| (t, p))
59            })
60            .collect::<OgeomResult<_>>()?;
61        for axis in AXES {
62            for sense in [1.0, -1.0] {
63                let score = |p: Point| p.to_vector().dot(axis) * sense;
64                let Some((best, _)) = samples
65                    .iter()
66                    .enumerate()
67                    .max_by(|a, b| score(a.1.1).total_cmp(&score(b.1.1)))
68                else {
69                    continue;
70                };
71                let lo = samples[best.saturating_sub(1)].0;
72                let hi = samples[(best + 1).min(N)].0;
73                points.push(samples[best].1);
74                if hi > lo {
75                    let found = solve::minimize(
76                        |t| at(t).map_or(f64::INFINITY, |p| -score(p)),
77                        lo,
78                        hi,
79                        solve::Criteria::default(),
80                    )?;
81                    points.push(at(found.value)?);
82                }
83            }
84        }
85    }
86
87    // Inside every face: the mesh's leading vertex for each side, moved on
88    // the exact surface as far as the side leads while it stays inside the
89    // face's own chart triangles.
90    for face in explore_unique(model, shape, ShapeType::Face)? {
91        let Some(NodeData::Face(data)) = model.node(&face).map(|n| n.data()) else {
92            continue;
93        };
94        let Some(surface) = model.geometry().surface(data.surface) else {
95            ogeom_bail!(Dangling, "surface is not in this model");
96        };
97        let placement = face.transform(model.datums())?;
98        let mesh = ogeom_mesh::triangulate_face(model, &face, Deflection::default(), tol)?;
99        if mesh.positions.is_empty() {
100            continue;
101        }
102        let chart: Vec<[Point2; 3]> = mesh
103            .triangles
104            .iter()
105            .map(|t| {
106                t.map(|i| {
107                    let (u, v) = mesh.parameters[i as usize];
108                    Point2::new(u, v)
109                })
110            })
111            .collect();
112        let inside = |p: Point2| chart.iter().any(|t| in_triangle(*t, p));
113        let span = chart.iter().flatten().fold(
114            (
115                f64::INFINITY,
116                f64::NEG_INFINITY,
117                f64::INFINITY,
118                f64::NEG_INFINITY,
119            ),
120            |b, p| (b.0.min(p.x), b.1.max(p.x), b.2.min(p.y), b.3.max(p.y)),
121        );
122        let lower = [span.0, span.2];
123        let upper = [span.1, span.3];
124        if !(upper[0] > lower[0] && upper[1] > lower[1]) {
125            continue;
126        }
127        for axis in AXES {
128            for sense in [1.0, -1.0] {
129                let score = |p: Point| p.to_vector().dot(axis) * sense;
130                let Some(best) = (0..mesh.positions.len())
131                    .max_by(|a, b| score(mesh.positions[*a]).total_cmp(&score(mesh.positions[*b])))
132                else {
133                    continue;
134                };
135                points.push(mesh.positions[best]);
136                let (u, v) = mesh.parameters[best];
137                let refined = solve_on_face(
138                    |uv: &[f64]| {
139                        let p = Point2::new(uv[0], uv[1]);
140                        if !inside(p) {
141                            return f64::INFINITY;
142                        }
143                        surface
144                            .point_at(uv[0], uv[1], tol)
145                            .map_or(f64::INFINITY, |q| -score(placement.apply(q)))
146                    },
147                    [u, v],
148                    lower,
149                    upper,
150                )?;
151                if let Some([u, v]) = refined {
152                    points.push(placement.apply(surface.point_at(u, v, tol)?));
153                }
154            }
155        }
156    }
157
158    if points.is_empty() {
159        ogeom_bail!(Construction, "the shape holds nothing to bound");
160    }
161    Ok(Aabb::of_points(&points))
162}
163
164/// A local descent from `start` inside the chart window; `None` where it
165/// found nothing better than the start.
166fn solve_on_face(
167    f: impl FnMut(&[f64]) -> f64,
168    start: [f64; 2],
169    lower: [f64; 2],
170    upper: [f64; 2],
171) -> OgeomResult<Option<[f64; 2]>> {
172    let step = ((upper[0] - lower[0]).min(upper[1] - lower[1]) * 1e-2).max(1e-9);
173    let found = ogeom_math::minimize_local(f, &start, &lower, &upper, step, 1e-15, 2000)?;
174    Ok(found
175        .value
176        .is_finite()
177        .then(|| [found.point[0], found.point[1]]))
178}
179
180fn in_triangle(t: [Point2; 3], p: Point2) -> bool {
181    let cross =
182        |a: Point2, b: Point2, c: Point2| (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
183    let (d1, d2, d3) = (
184        cross(t[0], t[1], p),
185        cross(t[1], t[2], p),
186        cross(t[2], t[0], p),
187    );
188    let eps = 1e-12;
189    let neg = d1 < -eps || d2 < -eps || d3 < -eps;
190    let pos = d1 > eps || d2 > eps || d3 > eps;
191    !(neg && pos)
192}
ogeom_algo/tight.rs

ogeom_algo/
tight.rs
