ogeom_algo/distance/
shape.rs

1//! Distances to a shape's own geometry.
2
3use std::sync::OnceLock;
4
5use ogeom_core::{OgeomError, OgeomResult, Tolerances, parallel};
6use ogeom_geom::Surface as _;
7use ogeom_math::{Aabb, Point, Point2, Vector};
8use ogeom_mesh::{Deflection, face_boundary_in_chart_met};
9use ogeom_topo::{Model, Orientation, SameKey, Shape, ShapeType, explore_unique};
10
11use super::{Closest, DistanceOptions, ShapeFeature, signed, unit, unit_or_zero};
12use crate::classify::{Containment, SolidBoundary, segment_distance};
13use crate::measure::{curve_bounds, project_on_curve, project_on_surface};
14use crate::proximity::{Element, PointDistance, RING_ANGULAR, RING_CHORD, Trim, trim_side};
15
16/// The signed distance from each point to `shape`'s faces, edges and
17/// vertices, measured on their own surfaces, curves and points, with the
18/// closest point, the feature it lies on and the normal there; in the
19/// points' order, `None` where nothing is within [`DistanceOptions::max`].
20///
21/// The candidates are [`PointDistance`]'s: each vertex, the foot on each
22/// edge, and the foot on each face's surface where it lands inside the
23/// face. A foot within the trimming rings' polylining band of the face's
24/// boundary is placed again against rings drawn a hundred times finer, and
25/// again, down to the confusion distance, for as long as it is nearer than
26/// every other candidate and the band cannot tell; one still in the band
27/// where the rings can be drawn no finer is left to the edges and vertices
28/// round it, as [`PointDistance::nearest`] leaves it. Equally near
29/// candidates go to the vertex before the edge before the face, and among
30/// those to the first the shape holds. A point that is not finite, or with
31/// a facing limit a normal that has no direction, answers `None`. See the
32/// [module](super) for how the sign is read.
33///
34/// A face's foot is its surface's nearest point to the query, as
35/// [`project_on_surface`] finds it: on a surface with more than one local
36/// nearest point, a nearer one inside the trim than the one found outside
37/// it is not a candidate.
38///
39/// # Errors
40///
41/// As [`to_mesh`](super::to_mesh) for the options and normals; as
42/// [`PointDistance::of`] for the shape, and as the projections and
43/// [`SolidBoundary::holds`] report.
44pub fn to_shape(
45    model: &Model,
46    shape: &Shape,
47    points: &[Point],
48    normals: Option<&[Vector]>,
49    options: &DistanceOptions,
50    tol: Tolerances,
51) -> OgeomResult<Vec<Option<Closest<ShapeFeature>>>> {
52    let facing = options.checked(points.len(), normals)?;
53    let target = Target::of(model, shape, options.signed, tol)?;
54    let answers = parallel::map_ordered(points, |i, &p| {
55        let wanted = match (facing, normals) {
56            (Some(cos), Some(normals)) => match unit(normals[i]) {
57                Some(m) => Some((m, cos)),
58                None => return Ok(None),
59            },
60            _ => None,
61        };
62        let faces = wanted.map(|(m, cos)| move |n: Vector| n.dot(m) >= cos);
63        let faces = faces.as_ref().map(|f| f as &dyn Fn(Vector) -> bool);
64        target.closest(model, p, faces, options, tol)
65    });
66    answers.into_iter().collect()
67}
68
69/// A face's trimming rings in its chart.
70type Rings = Vec<Vec<Point2>>;
71
72/// A shape's elements gathered once, with what each point asks of them.
73pub(crate) struct Target {
74    elements: PointDistance,
75    /// Each element's box, widened by the confusion distance.
76    bounds: Vec<Aabb>,
77    /// For an edge or a vertex, the faces of the shape that hold it, by
78    /// element index; empty for a face.
79    faces_of: Vec<Vec<usize>>,
80    /// For a face, its rings drawn at each of `chords`, as first needed;
81    /// `None` inside where they cannot be drawn within the chord.
82    finer: Vec<Vec<OnceLock<Option<Rings>>>>,
83    /// The chords rings are drawn again at, finest last.
84    chords: Vec<f64>,
85    /// The shape's boundary prepared for classifying, where it is closed
86    /// and the distances are signed.
87    solid: Option<SolidBoundary>,
88}
89
90/// What a nearest point passes over: the elements held by `faces` alone
91/// (a face of them, or an edge or a vertex every face holding which is one
92/// of them) whose nearest point lies within the confusion distance of the
93/// segment from `from` to the point asked about.
94pub(crate) struct Passed<'a> {
95    pub(crate) faces: &'a [Shape],
96    pub(crate) from: Point,
97}
98
99/// One candidate: its unsigned distance, its point, its element, and the
100/// normal its side is read from.
101#[derive(Debug, Clone, Copy)]
102pub(crate) struct Found {
103    pub(crate) distance: f64,
104    pub(crate) point: Point,
105    element: usize,
106    normal: Vector,
107}
108
109impl Target {
110    /// The shape's elements gathered, with its boundary prepared for
111    /// classifying where `sign` asks for signed distances.
112    pub(crate) fn of(
113        model: &Model,
114        shape: &Shape,
115        sign: bool,
116        tol: Tolerances,
117    ) -> OgeomResult<Self> {
118        let elements = PointDistance::of(model, shape, tol)?;
119        let list = elements.elements();
120        let mut places: ogeom_core::FastMap<SameKey, usize> = ogeom_core::FastMap::default();
121        let mut bounds = Vec::with_capacity(list.len());
122        let reach = tol.confusion().max(model.tolerances().confusion());
123        for (i, element) in list.iter().enumerate() {
124            let (held, bound) = match element {
125                Element::Vertex(v, p) => (v, Aabb::of_point(*p)),
126                Element::Edge(e, curve) => (e, curve_bounds(curve, tol)?),
127                Element::Face(f, _) => (f, crate::tight::face_bounds(model, f)?),
128            };
129            places.insert(SameKey(held.clone()), i);
130            bounds.push(bound.expanded(reach));
131        }
132        let mut faces_of = vec![Vec::new(); list.len()];
133        let mut finer = Vec::with_capacity(list.len());
134        let mut chords = Vec::new();
135        let mut chord = RING_CHORD;
136        while chord > tol.confusion() {
137            chord = (chord / 100.0).max(tol.confusion());
138            chords.push(chord);
139        }
140        for (i, element) in list.iter().enumerate() {
141            let Element::Face(face, _) = element else {
142                finer.push(Vec::new());
143                continue;
144            };
145            finer.push(chords.iter().map(|_| OnceLock::new()).collect());
146            for kind in [ShapeType::Edge, ShapeType::Vertex] {
147                for held in explore_unique(model, face, kind)? {
148                    if let Some(&at) = places.get(&SameKey(held)) {
149                        faces_of[at].push(i);
150                    }
151                }
152            }
153        }
154        let solid = if sign {
155            let chord = crate::classify::exact_ring_chord(tol);
156            match SolidBoundary::of(model, shape, chord, tol) {
157                Ok(solid) => Some(solid),
158                // An open boundary has no inside to classify against.
159                Err(OgeomError::Construction(_)) => None,
160                Err(other) => return Err(other),
161            }
162        } else {
163            None
164        };
165        Ok(Self {
166            elements,
167            bounds,
168            faces_of,
169            finer,
170            chords,
171            solid,
172        })
173    }
174
175    /// The oriented unit normal of face element `face` at `(u, v)` on its
176    /// world surface, or `None` where the surface has none.
177    fn normal(&self, face: usize, (u, v): (f64, f64), tol: Tolerances) -> Option<Vector> {
178        let Element::Face(shape, prepared) = &self.elements.elements()[face] else {
179            return None;
180        };
181        let n = Vector::from(prepared.world.normal_at(u, v, tol).ok()?);
182        Some(if shape.orientation() == Orientation::Reversed {
183            -n
184        } else {
185            n
186        })
187    }
188
189    /// The unit normals at `q`, on the faces that hold edge or vertex
190    /// element `element`, of those that have one there.
191    fn normals_at(&self, element: usize, q: Point, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
192        let mut out = Vec::with_capacity(self.faces_of[element].len());
193        for &face in &self.faces_of[element] {
194            let Element::Face(_, prepared) = &self.elements.elements()[face] else {
195                continue;
196            };
197            let foot = project_on_surface(&prepared.world, q, 32, tol)?;
198            out.extend(self.normal(face, foot.parameters, tol));
199        }
200        Ok(out)
201    }
202
203    fn closest(
204        &self,
205        model: &Model,
206        p: Point,
207        faces: Option<&dyn Fn(Vector) -> bool>,
208        options: &DistanceOptions,
209        tol: Tolerances,
210    ) -> OgeomResult<Option<Closest<ShapeFeature>>> {
211        if !(p.x.is_finite() && p.y.is_finite() && p.z.is_finite()) {
212            return Ok(None);
213        }
214        let Some(found) = self.nearest(model, p, faces, None, options.max, tol)? else {
215            return Ok(None);
216        };
217        let element = &self.elements.elements()[found.element];
218        let mut normal = found.normal;
219        let feature = match element {
220            Element::Vertex(v, _) => ShapeFeature::Vertex(v.clone()),
221            Element::Edge(e, _) => ShapeFeature::Edge(e.clone()),
222            Element::Face(f, _) => ShapeFeature::Face(f.clone()),
223        };
224        let mut distance = signed(found.distance, p, found.point, normal, options.signed);
225        if options.signed
226            && matches!(feature, ShapeFeature::Vertex(_))
227            && let Some(solid) = &self.solid
228        {
229            match solid.holds(model, p, tol)? {
230                Containment::In => distance = -found.distance,
231                Containment::Out => distance = found.distance,
232                Containment::On => {}
233            }
234        }
235        normal = unit_or_zero(normal);
236        Ok(Some(Closest {
237            distance,
238            point: found.point,
239            feature,
240            normal,
241        }))
242    }
243
244    /// The nearest candidate within `max` whose face's unit normal there
245    /// `faces` accepts, where there is a filter (for an edge or a vertex,
246    /// the normal of any face holding it): elements taken nearest box first
247    /// and passed over once their box is farther than the best so far,
248    /// which changes only the time. The elements `passed` names are passed
249    /// over.
250    pub(crate) fn nearest(
251        &self,
252        model: &Model,
253        p: Point,
254        faces: Option<&dyn Fn(Vector) -> bool>,
255        passed: Option<&Passed<'_>>,
256        max: f64,
257        tol: Tolerances,
258    ) -> OgeomResult<Option<Found>> {
259        let list = self.elements.elements();
260        let mut order: Vec<(f64, usize)> = self
261            .bounds
262            .iter()
263            .enumerate()
264            .map(|(i, b)| (b.distance_to(p), i))
265            .filter(|&(gap, _)| gap <= max)
266            .collect();
267        order.sort_by(|a, b| a.0.total_cmp(&b.0).then(a.1.cmp(&b.1)));
268        let mut best: Option<Found> = None;
269        let better = |best: &Option<Found>, distance: f64, element: usize| match best {
270            None => distance <= max,
271            Some(held) => {
272                distance < held.distance || (distance == held.distance && element < held.element)
273            }
274        };
275        let kept = |normals: &[Vector]| faces.is_none_or(|faces| normals.iter().any(|&n| faces(n)));
276        let skipped = |element: usize, q: Point| {
277            passed.is_some_and(|passed| {
278                segment_distance(q, passed.from, p) <= tol.confusion()
279                    && self.held_by(element, passed.faces)
280            })
281        };
282        let sum = |normals: Vec<Vector>| normals.into_iter().fold(Vector::ZERO, |a, n| a + n);
283        let mut banded: Vec<Found> = Vec::new();
284        for (gap, i) in order {
285            if best.is_some_and(|held| gap > held.distance) {
286                break;
287            }
288            match &list[i] {
289                Element::Vertex(_, q) => {
290                    let distance = p.distance(*q);
291                    if !better(&best, distance, i) || skipped(i, *q) {
292                        continue;
293                    }
294                    let normals = self.normals_at(i, *q, tol)?;
295                    if kept(&normals) {
296                        best = Some(Found {
297                            distance,
298                            point: *q,
299                            element: i,
300                            normal: sum(normals),
301                        });
302                    }
303                }
304                Element::Edge(_, curve) => {
305                    let foot = project_on_curve(curve, p, 64, tol)?;
306                    if !better(&best, foot.distance, i) || skipped(i, foot.point) {
307                        continue;
308                    }
309                    let normals = self.normals_at(i, foot.point, tol)?;
310                    if kept(&normals) {
311                        best = Some(Found {
312                            distance: foot.distance,
313                            point: foot.point,
314                            element: i,
315                            normal: sum(normals),
316                        });
317                    }
318                }
319                Element::Face(_, prepared) => {
320                    let foot = project_on_surface(&prepared.world, p, 32, tol)?;
321                    if !better(&best, foot.distance, i) || skipped(i, foot.point) {
322                        continue;
323                    }
324                    // A foot where the surface has no normal (a pole, an
325                    // apex) is a point of the edges and vertices there.
326                    let Some(normal) = self.normal(i, foot.parameters, tol) else {
327                        continue;
328                    };
329                    if !faces.is_none_or(|faces| faces(normal)) {
330                        continue;
331                    }
332                    let found = Found {
333                        distance: foot.distance,
334                        point: foot.point,
335                        element: i,
336                        normal,
337                    };
338                    match trim_side(prepared, None, RING_CHORD, foot.point, tol)? {
339                        Trim::Inside => best = Some(found),
340                        Trim::Outside => {}
341                        Trim::Band => banded.push(found),
342                    }
343                }
344            }
345        }
346        banded.sort_by(|a, b| {
347            a.distance
348                .total_cmp(&b.distance)
349                .then(a.element.cmp(&b.element))
350        });
351        for found in banded {
352            if !better(&best, found.distance, found.element) {
353                break;
354            }
355            if self.inside_finer(model, found.element, found.point, tol)? {
356                best = Some(found);
357                break;
358            }
359        }
360        Ok(best)
361    }
362
363    /// Whether element `element` is held by `faces` alone: is one of them,
364    /// or is an edge or a vertex every face holding which is.
365    fn held_by(&self, element: usize, faces: &[Shape]) -> bool {
366        let list = self.elements.elements();
367        let among = |face: usize| match &list[face] {
368            Element::Face(shape, _) => faces.iter().any(|own| own.is_same(shape)),
369            _ => false,
370        };
371        match &list[element] {
372            Element::Face(..) => among(element),
373            _ => {
374                let holding = &self.faces_of[element];
375                !holding.is_empty() && holding.iter().all(|&face| among(face))
376            }
377        }
378    }
379
380    /// Whether a foot the coarse rings put in their band is inside its
381    /// face, placed against rings drawn finer until it leaves the band;
382    /// `false` where it is outside, or still in the band where the rings
383    /// can be drawn no finer.
384    fn inside_finer(
385        &self,
386        model: &Model,
387        face: usize,
388        q: Point,
389        tol: Tolerances,
390    ) -> OgeomResult<bool> {
391        let Element::Face(shape, prepared) = &self.elements.elements()[face] else {
392            return Ok(false);
393        };
394        for (level, &chord) in self.chords.iter().enumerate() {
395            let rings = self.finer[face][level].get_or_init(|| {
396                // A curve's segments at a chord grow as the chord's inverse
397                // square root, and so does the cap: the coarse rings' cap
398                // scaled by as much.
399                let coarse = Deflection::default().max_segments;
400                #[allow(
401                    clippy::cast_possible_truncation,
402                    clippy::cast_sign_loss,
403                    clippy::cast_precision_loss
404                )]
405                let cap = (coarse as f64 * (RING_CHORD / chord).sqrt()).ceil() as usize;
406                let deflection = Deflection {
407                    chord,
408                    angular: RING_ANGULAR,
409                    max_segments: cap,
410                    ..Deflection::default()
411                };
412                // A boundary that will not draw finer leaves the foot in
413                // the band, where the coarse rings left it.
414                face_boundary_in_chart_met(model, shape, deflection, tol)
415                    .ok()
416                    .flatten()
417            });
418            let Some(rings) = rings else {
419                return Ok(false);
420            };
421            match trim_side(prepared, Some(rings), chord, q, tol)? {
422                Trim::Inside => return Ok(true),
423                Trim::Outside => return Ok(false),
424                Trim::Band => {}
425            }
426        }
427        Ok(false)
428    }
429}
ogeom_algo/distance/shape.rs

ogeom_algo/distance/
shape.rs
