axiolid_spatial/
barycentric.rs

1//! Barycentric and mean-value coordinates, for interpolating values given at
2//! the corners of a triangle, a tetrahedron or a polygon (#143).
3//!
4//! Every function returns one weight per corner. The weights sum to one and
5//! reproduce the point: `sum(w_i * corner_i) == point`, up to rounding. A
6//! point on a corner gets weight exactly one there and zero elsewhere, so
7//! corner values are reproduced exactly.
8//!
9//! Shapes thinner than the linear tolerance are refused rather than given
10//! weights that are mostly rounding noise; a triangle's or tetrahedron's
11//! thinness is its least altitude. Weights outside `[0, 1]` are not an
12//! error: they place the point outside the triangle or tetrahedron, which is
13//! how callers extrapolate or test containment.
14
15use axiolid_core::{Point2, Point3, Polygon2, Scalar, Tolerance, Triangle2, Triangle3};
16
17/// Why coordinates could not be computed.
18#[derive(Debug, Clone, Copy, PartialEq)]
19#[non_exhaustive]
20pub enum BarycentricError {
21    /// A corner or the query point is not finite.
22    NonFinite,
23    /// The triangle, tetrahedron or polygon is thinner than the linear
24    /// tolerance: its least altitude (for a polygon, twice its area over its
25    /// perimeter) is `thickness`.
26    Degenerate {
27        /// The shape's thickness, in length units.
28        thickness: Scalar,
29    },
30    /// A polygon has fewer than three vertices.
31    TooFewVertices {
32        /// The number of vertices given.
33        count: usize,
34    },
35    /// Polygon edge `index`, from vertex `index` to the next, is no longer
36    /// than the linear tolerance.
37    ShortEdge {
38        /// The edge's index.
39        index: usize,
40    },
41    /// Polygon edges `first` and `second` cross, touch or fold back onto
42    /// each other within the linear tolerance, so the polygon is not simple.
43    SelfIntersecting {
44        /// The lower edge index.
45        first: usize,
46        /// The higher edge index.
47        second: usize,
48    },
49    /// The mean-value weights cancel at this point, which can happen only
50    /// outside a non-convex polygon: the coordinates are undefined there.
51    Undefined,
52}
53
54impl core::fmt::Display for BarycentricError {
55    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
56        match self {
57            Self::NonFinite => f.write_str("a corner or the query point is not finite"),
58            Self::Degenerate { thickness } => write!(
59                f,
60                "shape is {thickness} thick, not thicker than the linear tolerance"
61            ),
62            Self::TooFewVertices { count } => {
63                write!(f, "polygon has {count} vertices, need at least 3")
64            }
65            Self::ShortEdge { index } => {
66                write!(f, "polygon edge {index} is not longer than the tolerance")
67            }
68            Self::SelfIntersecting { first, second } => {
69                write!(f, "polygon edges {first} and {second} meet")
70            }
71            Self::Undefined => {
72                f.write_str("mean-value weights cancel here (outside a non-convex polygon)")
73            }
74        }
75    }
76}
77
78impl core::error::Error for BarycentricError {}
79
80/// Twice the signed area of `x, y, z`, positive counter-clockwise.
81fn area2(x: Point2, y: Point2, z: Point2) -> Scalar {
82    (y - x).perp_dot(z - x)
83}
84
85/// Six times the signed volume of `x, y, z, w`.
86fn volume6(x: Point3, y: Point3, z: Point3, w: Point3) -> Scalar {
87    (y - x).dot((z - x).cross(w - x))
88}
89
90/// Barycentric coordinates of `point` in a 2D triangle, as weights of
91/// `a`, `b` and `c`.
92///
93/// Each weight is the signed area of the sub-triangle opposite its corner
94/// over the triangle's area, so the result does not depend on the
95/// triangle's winding.
96///
97/// # Errors
98///
99/// [`BarycentricError::NonFinite`], or [`BarycentricError::Degenerate`] when
100/// the triangle's least altitude is not above `tolerance.linear()`.
101pub fn triangle_barycentric2(
102    triangle: &Triangle2,
103    point: Point2,
104    tolerance: Tolerance,
105) -> Result<[Scalar; 3], BarycentricError> {
106    let Triangle2 { a, b, c } = *triangle;
107    if ![a, b, c, point].iter().all(|p| p.is_finite()) {
108        return Err(BarycentricError::NonFinite);
109    }
110    let total = area2(a, b, c);
111    let longest = (b - a).length().max((c - b).length()).max((a - c).length());
112    thick_enough(total.abs(), longest, tolerance)?;
113    Ok([
114        area2(point, b, c) / total,
115        area2(a, point, c) / total,
116        area2(a, b, point) / total,
117    ])
118}
119
120/// Barycentric coordinates of `point` in a 3D triangle, as weights of
121/// `a`, `b` and `c`.
122///
123/// A point off the triangle's plane gets the coordinates of its orthogonal
124/// projection onto that plane: the weights reproduce the projection, not the
125/// point.
126///
127/// # Errors
128///
129/// As [`triangle_barycentric2`].
130pub fn triangle_barycentric3(
131    triangle: &Triangle3,
132    point: Point3,
133    tolerance: Tolerance,
134) -> Result<[Scalar; 3], BarycentricError> {
135    let Triangle3 { a, b, c } = *triangle;
136    if ![a, b, c, point].iter().all(|p| p.is_finite()) {
137        return Err(BarycentricError::NonFinite);
138    }
139    let normal = (b - a).cross(c - a);
140    let longest = (b - a).length().max((c - b).length()).max((a - c).length());
141    thick_enough(normal.length(), longest, tolerance)?;
142    let total = normal.dot(normal);
143    Ok([
144        normal.dot((b - point).cross(c - point)) / total,
145        normal.dot((c - point).cross(a - point)) / total,
146        normal.dot((a - point).cross(b - point)) / total,
147    ])
148}
149
150/// Barycentric coordinates of `point` in a tetrahedron, as weights of its
151/// four corners in order.
152///
153/// Each weight is the signed volume of the sub-tetrahedron opposite its
154/// corner over the tetrahedron's volume, so the corner order's handedness
155/// does not matter.
156///
157/// # Errors
158///
159/// [`BarycentricError::NonFinite`], or [`BarycentricError::Degenerate`] when
160/// the tetrahedron's least altitude is not above `tolerance.linear()`.
161pub fn tetrahedron_barycentric(
162    corners: [Point3; 4],
163    point: Point3,
164    tolerance: Tolerance,
165) -> Result<[Scalar; 4], BarycentricError> {
166    let [a, b, c, d] = corners;
167    if ![a, b, c, d, point].iter().all(|p| p.is_finite()) {
168        return Err(BarycentricError::NonFinite);
169    }
170    let total = volume6(a, b, c, d);
171    // The least altitude is the volume over the largest face.
172    let largest_face = [(b, c, d), (a, c, d), (a, b, d), (a, b, c)]
173        .iter()
174        .map(|&(x, y, z)| (y - x).cross(z - x).length())
175        .fold(0.0, Scalar::max);
176    thick_enough(total.abs(), largest_face, tolerance)?;
177    Ok([
178        volume6(point, b, c, d) / total,
179        volume6(a, point, c, d) / total,
180        volume6(a, b, point, d) / total,
181        volume6(a, b, c, point) / total,
182    ])
183}
184
185/// Refuse a shape whose `measure / base` (an altitude) is not above the
186/// linear tolerance.
187fn thick_enough(
188    measure: Scalar,
189    base: Scalar,
190    tolerance: Tolerance,
191) -> Result<(), BarycentricError> {
192    let thickness = if base > 0.0 { measure / base } else { 0.0 };
193    if thickness > tolerance.linear() {
194        Ok(())
195    } else {
196        Err(BarycentricError::Degenerate { thickness })
197    }
198}
199
200/// Mean-value coordinates of `point` in a simple polygon, one weight per
201/// vertex (Floater 2003, in the form for arbitrary polygons of Hormann and
202/// Floater 2006).
203///
204/// Inside a simple polygon, convex or not, the weights are positive-sum,
205/// smooth, and reproduce the point. On the boundary, within the linear
206/// tolerance, they are the boundary's own linear interpolation: one at a
207/// vertex, or the two endpoint weights of an edge. Outside the polygon they
208/// are defined wherever the raw weights do not cancel, which is everywhere
209/// outside a convex polygon. The polygon may wind either way.
210///
211/// The polygon is checked for simplicity pairwise, `O(n^2)` in its vertex
212/// count, since mean-value coordinates of a crossing polygon mean nothing.
213///
214/// # Errors
215///
216/// [`BarycentricError::NonFinite`], [`BarycentricError::TooFewVertices`],
217/// [`BarycentricError::ShortEdge`], [`BarycentricError::SelfIntersecting`],
218/// [`BarycentricError::Degenerate`] for a polygon with no area to speak of,
219/// and [`BarycentricError::Undefined`] where the weights cancel.
220pub fn mean_value_coordinates2(
221    polygon: &Polygon2,
222    point: Point2,
223    tolerance: Tolerance,
224) -> Result<Vec<Scalar>, BarycentricError> {
225    let v = &polygon.vertices;
226    let n = v.len();
227    if !point.is_finite() || !v.iter().all(|p| p.is_finite()) {
228        return Err(BarycentricError::NonFinite);
229    }
230    if n < 3 {
231        return Err(BarycentricError::TooFewVertices { count: n });
232    }
233    check_simple(v, tolerance)?;
234    let perimeter: Scalar = (0..n).map(|i| (v[(i + 1) % n] - v[i]).length()).sum();
235    thick_enough(2.0 * polygon.signed_area().abs(), perimeter, tolerance)?;
236
237    let linear = tolerance.linear();
238    let mut weights = vec![0.0; n];
239    // On a vertex: that vertex alone.
240    let s: Vec<Point2> = v.iter().map(|&q| q - point).collect();
241    let r: Vec<Scalar> = s.iter().map(|q| q.length()).collect();
242    if let Some(i) = (0..n)
243        .filter(|&i| r[i] <= linear)
244        .min_by(|&i, &j| r[i].total_cmp(&r[j]))
245    {
246        weights[i] = 1.0;
247        return Ok(weights);
248    }
249    // On an edge: its linear interpolation.
250    for i in 0..n {
251        let j = (i + 1) % n;
252        let edge = v[j] - v[i];
253        let length = edge.length();
254        let t = (point - v[i]).dot(edge) / (length * length);
255        if (0.0..=1.0).contains(&t) && s[i].perp_dot(s[j]).abs() / length <= linear {
256            weights[i] = 1.0 - t;
257            weights[j] = t;
258            return Ok(weights);
259        }
260    }
261    // tan(alpha_i / 2) for the signed angle alpha_i the edge from vertex i
262    // to vertex i + 1 subtends at the point: A / (r_i r_j + D), which stays
263    // finite because the point is on no edge.
264    let half_tangent: Vec<Scalar> = (0..n)
265        .map(|i| {
266            let j = (i + 1) % n;
267            s[i].perp_dot(s[j]) / (r[i] * r[j] + s[i].dot(s[j]))
268        })
269        .collect();
270    let mut sum = 0.0;
271    let mut magnitude = 0.0;
272    for i in 0..n {
273        let w = (half_tangent[(i + n - 1) % n] + half_tangent[i]) / r[i];
274        weights[i] = w;
275        sum += w;
276        magnitude += w.abs();
277    }
278    if !sum.is_finite() || sum.abs() <= Scalar::EPSILON * n as Scalar * magnitude {
279        return Err(BarycentricError::Undefined);
280    }
281    for w in &mut weights {
282        *w /= sum;
283    }
284    Ok(weights)
285}
286
287/// Refuse short edges, and non-adjacent edges that come within the linear
288/// tolerance of each other.
289///
290/// Adjacent edges need no test of their own. Two that fold back onto each
291/// other leave the far end of the shorter on the longer, and that end also
292/// belongs to an edge not adjacent to the longer, which this test sees; in a
293/// triangle, where every pair is adjacent, a fold leaves no area and the
294/// thickness test refuses it.
295fn check_simple(v: &[Point2], tolerance: Tolerance) -> Result<(), BarycentricError> {
296    let n = v.len();
297    let linear = tolerance.linear();
298    let edge = |i: usize| (v[i], v[(i + 1) % n]);
299    for i in 0..n {
300        let (p, q) = edge(i);
301        if (q - p).length() <= linear {
302            return Err(BarycentricError::ShortEdge { index: i });
303        }
304    }
305    for i in 0..n {
306        // Skip the edge itself, its successor, and (from edge 0) the last
307        // edge, which precedes it.
308        let last = if i == 0 { n - 1 } else { n };
309        for j in i + 2..last {
310            let (p, q) = edge(i);
311            let (r, s) = edge(j);
312            if segment_distance(p, q, r, s) <= linear {
313                return Err(BarycentricError::SelfIntersecting {
314                    first: i,
315                    second: j,
316                });
317            }
318        }
319    }
320    Ok(())
321}
322
323/// Distance between segments `pq` and `rs`: zero when they cross.
324fn segment_distance(p: Point2, q: Point2, r: Point2, s: Point2) -> Scalar {
325    let o1 = area2(p, q, r);
326    let o2 = area2(p, q, s);
327    let o3 = area2(r, s, p);
328    let o4 = area2(r, s, q);
329    if o1 * o2 < 0.0 && o3 * o4 < 0.0 {
330        return 0.0;
331    }
332    point_segment(r, p, q)
333        .min(point_segment(s, p, q))
334        .min(point_segment(p, r, s))
335        .min(point_segment(q, r, s))
336}
337
338fn point_segment(x: Point2, p: Point2, q: Point2) -> Scalar {
339    let d = q - p;
340    let t = ((x - p).dot(d) / d.dot(d)).clamp(0.0, 1.0);
341    (p + d * t - x).length()
342}
axiolid_spatial/barycentric.rs

axiolid_spatial/
barycentric.rs
