geop_ops_rasterize/
polygon_triangulate.rs

1//! 2-D ear-clipping triangulation of a simple polygon with holes.
2//!
3//! Used to rasterize trimmed faces (faces with more than one edge_loop): such faces are
4//! planar/affine in this codebase, so triangulating the `(u, v)`
5//! outer-loop-minus-holes polygon and mapping each triangle through
6//! `surface.eval_at` gives an exact 3-D triangulation.
7
8use geop_core_math::{geop_error::GeopResult, scalars::Scalar, vector::Vector2};
9
10pub use geop_core_math::polygon::polygon_signed_area;
11
12/// Twice the signed area of `poly` (shoelace formula): positive for CCW,
13/// negative for CW.
14fn signed_area2<S: Scalar>(poly: &[Vector2<S>]) -> S {
15    let n = poly.len();
16    let mut sum = S::ZERO;
17    for i in 0..n {
18        let j = (i + 1) % n;
19        sum = sum.add(poly[i][0].mul(poly[j][1]).sub(poly[j][0].mul(poly[i][1])));
20    }
21    sum
22}
23
24/// 2-D cross product (perp-dot) of `a` and `b`.
25fn cross2<S: Scalar>(a: Vector2<S>, b: Vector2<S>) -> S {
26    a[0].mul(b[1]).sub(a[1].mul(b[0]))
27}
28
29fn to_f64_2<S: Scalar>(p: Vector2<S>) -> (f64, f64) {
30    (p[0].to_f64(), p[1].to_f64())
31}
32
33/// True if `p` lies inside (or on the edge_loop of) triangle `(a, b, c)`,
34/// using same-sign barycentric tests. Ambiguous (interval-straddling) sign
35/// comparisons are treated as "inside": this is conservative for ear-clipping
36/// (an ambiguous result rejects the candidate ear rather than risking an
37/// invalid triangulation).
38fn point_in_triangle<S: Scalar>(
39    p: Vector2<S>,
40    a: Vector2<S>,
41    b: Vector2<S>,
42    c: Vector2<S>,
43) -> bool {
44    let d1 = cross2(b.sub(&a), p.sub(&a));
45    let d2 = cross2(c.sub(&b), p.sub(&b));
46    let d3 = cross2(a.sub(&c), p.sub(&c));
47
48    let all_nonneg = !d1.definitely_less(S::ZERO)
49        && !d2.definitely_less(S::ZERO)
50        && !d3.definitely_less(S::ZERO);
51    let all_nonpos = !d1.definitely_greater(S::ZERO)
52        && !d2.definitely_greater(S::ZERO)
53        && !d3.definitely_greater(S::ZERO);
54    all_nonneg || all_nonpos
55}
56
57/// Ear-clipping triangulation of a simple polygon (no holes), returned as
58/// index triples into `poly`. Robust to interval scalars: convexity and
59/// point-in-triangle tests use three-valued comparisons, treating ambiguous
60/// results conservatively (reject the candidate ear); a fallback guarantees
61/// termination on degenerate input.
62fn ear_clip<S: Scalar>(poly: &[Vector2<S>]) -> Vec<(usize, usize, usize)> {
63    let n = poly.len();
64    if n < 3 {
65        return Vec::new();
66    }
67
68    let total = signed_area2(poly);
69    // Orientation used to decide which turn-sign counts as "convex". Default
70    // to CCW (positive) if the total area is ambiguous (near-zero).
71    let orient_positive = !total.definitely_less(S::ZERO);
72
73    let mut indices: Vec<usize> = (0..n).collect();
74    let mut triangles = Vec::with_capacity(n.saturating_sub(2));
75
76    let max_iters = n * n + 8;
77    let mut iters = 0;
78    while indices.len() > 3 && iters < max_iters {
79        iters += 1;
80        let m = indices.len();
81        let mut clipped = false;
82        for k in 0..m {
83            let ip = indices[(k + m - 1) % m];
84            let ic = indices[k];
85            let inext = indices[(k + 1) % m];
86
87            let turn = cross2(poly[ic].sub(&poly[ip]), poly[inext].sub(&poly[ic]));
88            let is_convex = if orient_positive {
89                turn.definitely_greater(S::ZERO)
90            } else {
91                turn.definitely_less(S::ZERO)
92            };
93            if !is_convex {
94                continue;
95            }
96
97            // A vertex that coincides with one of the triangle's own corners
98            // (e.g. the duplicated bridge endpoints introduced by
99            // `merge_hole_into`'s zero-width slit) lies on the triangle's
100            // edge_loop by construction, not in its interior — don't let that
101            // reject an otherwise-valid ear.
102            let contains_other = indices.iter().any(|&iq| {
103                iq != ip
104                    && iq != ic
105                    && iq != inext
106                    && !points_equal(poly[iq], poly[ip])
107                    && !points_equal(poly[iq], poly[ic])
108                    && !points_equal(poly[iq], poly[inext])
109                    && point_in_triangle(poly[iq], poly[ip], poly[ic], poly[inext])
110            });
111            if contains_other {
112                continue;
113            }
114
115            triangles.push((ip, ic, inext));
116            indices.remove(k);
117            clipped = true;
118            break;
119        }
120        if !clipped {
121            // Degenerate/ambiguous configuration: clip the first vertex
122            // regardless, to guarantee termination. `TriangleFace::try_new`
123            // (in the caller) will discard genuinely-degenerate triangles.
124            let m = indices.len();
125            let ip = indices[m - 1];
126            let ic = indices[0];
127            let inext = indices[1];
128            triangles.push((ip, ic, inext));
129            indices.remove(0);
130        }
131    }
132    if indices.len() == 3 {
133        triangles.push((indices[0], indices[1], indices[2]));
134    }
135    triangles
136}
137
138/// Squared Euclidean distance between `a` and `b`, via `to_f64()`. Used only
139/// to heuristically order bridge candidates — does not affect correctness.
140fn dist2_f64<S: Scalar>(a: Vector2<S>, b: Vector2<S>) -> f64 {
141    let (ax, ay) = to_f64_2(a);
142    let (bx, by) = to_f64_2(b);
143    let dx = ax - bx;
144    let dy = ay - by;
145    dx * dx + dy * dy
146}
147
148fn points_equal<S: Scalar>(a: Vector2<S>, b: Vector2<S>) -> bool {
149    a.could_be_equal(&b)
150}
151
152/// True if segments `(p1,p2)` and `(p3,p4)` properly cross (each segment's
153/// endpoints lie strictly on opposite sides of the other segment's line).
154/// Ambiguous (collinear/touching) configurations are treated as
155/// non-intersecting — shared endpoints are handled separately by the caller.
156fn segments_properly_intersect<S: Scalar>(
157    p1: Vector2<S>,
158    p2: Vector2<S>,
159    p3: Vector2<S>,
160    p4: Vector2<S>,
161) -> bool {
162    let d1 = cross2(p2.sub(&p1), p3.sub(&p1));
163    let d2 = cross2(p2.sub(&p1), p4.sub(&p1));
164    let d3 = cross2(p4.sub(&p3), p1.sub(&p3));
165    let d4 = cross2(p4.sub(&p3), p2.sub(&p3));
166
167    let opposite = |x: S, y: S| {
168        (x.definitely_greater(S::ZERO) && y.definitely_less(S::ZERO))
169            || (x.definitely_less(S::ZERO) && y.definitely_greater(S::ZERO))
170    };
171
172    opposite(d1, d2) && opposite(d3, d4)
173}
174
175/// Even-odd point-in-polygon test via `to_f64()` ray casting. Used only to
176/// pick a good bridge target when merging a hole into the outer polygon — not
177/// correctness-critical (the final ear-clipping result is validated with
178/// exact comparisons).
179fn point_in_polygon<S: Scalar>(p: Vector2<S>, poly: &[Vector2<S>]) -> bool {
180    let n = poly.len();
181    let (px, py) = to_f64_2(p);
182    let mut inside = false;
183    for i in 0..n {
184        let (ax, ay) = to_f64_2(poly[i]);
185        let (bx, by) = to_f64_2(poly[(i + 1) % n]);
186        if (ay > py) != (by > py) {
187            let x_intersect = ax + (py - ay) / (by - ay) * (bx - ax);
188            if px < x_intersect {
189                inside = !inside;
190            }
191        }
192    }
193    inside
194}
195
196/// True if the bridge segment `h -> m` (a candidate connection from a hole
197/// vertex to an outer-polygon vertex) doesn't cross `poly`'s edge_loop and its
198/// midpoint stays inside `poly`.
199fn is_bridge_visible<S: Scalar>(poly: &[Vector2<S>], h: Vector2<S>, m: Vector2<S>) -> bool {
200    if points_equal(h, m) {
201        return false;
202    }
203    let n = poly.len();
204    for i in 0..n {
205        let a = poly[i];
206        let b = poly[(i + 1) % n];
207        if points_equal(a, h) || points_equal(b, h) || points_equal(a, m) || points_equal(b, m) {
208            continue;
209        }
210        if segments_properly_intersect(h, m, a, b) {
211            return false;
212        }
213    }
214    let two = S::TWO;
215    let mid = Vector2::from_array([
216        h[0].add(m[0]).div(two).unwrap_or(h[0]),
217        h[1].add(m[1]).div(two).unwrap_or(h[1]),
218    ]);
219    point_in_polygon(mid, poly)
220}
221
222/// Merge `hole` (already oriented CW) into `merged` (already oriented CCW) by
223/// finding a visible bridge and splicing the hole's vertices in, producing a
224/// single simple polygon with a zero-width slit.
225fn merge_hole_into<S: Scalar>(merged: &mut Vec<Vector2<S>>, hole: &[Vector2<S>]) {
226    // Bridge start: the hole vertex with max x (tie-break max y). This choice
227    // only affects which bridge is found, not correctness.
228    let hi = (1..hole.len()).fold(0usize, |best, i| {
229        let (bx, by) = to_f64_2(hole[best]);
230        let (ix, iy) = to_f64_2(hole[i]);
231        if ix > bx || (ix == bx && iy > by) {
232            i
233        } else {
234            best
235        }
236    });
237    let h = hole[hi];
238
239    let mut candidates: Vec<usize> = (0..merged.len()).collect();
240    candidates.sort_by(|&a, &b| {
241        dist2_f64(h, merged[a])
242            .partial_cmp(&dist2_f64(h, merged[b]))
243            .unwrap()
244    });
245
246    let mi = candidates
247        .into_iter()
248        .find(|&i| is_bridge_visible(merged, h, merged[i]))
249        .unwrap_or(0);
250    let m = merged[mi];
251
252    let mut new_merged = Vec::with_capacity(merged.len() + hole.len() + 2);
253    new_merged.extend_from_slice(&merged[0..=mi]);
254    for k in 0..hole.len() {
255        new_merged.push(hole[(hi + k) % hole.len()]);
256    }
257    new_merged.push(h);
258    new_merged.push(m);
259    new_merged.extend_from_slice(&merged[mi + 1..]);
260
261    *merged = new_merged;
262}
263
264/// Merge `outer` (a simple polygon) and `holes` (simple polygons strictly
265/// inside `outer` and disjoint from each other) into one simple polygon —
266/// each hole spliced in via a zero-width slit bridge (see
267/// [`merge_hole_into`]) — oriented CCW with every hole CW. The result
268/// traces exactly the trimmed region's boundary, so any triangulator (or,
269/// in `adaptive`, a constrained-Delaunay triangulator) that respects it as
270/// a closed ring automatically excludes the holes without further
271/// bookkeeping.
272pub(crate) fn merge_outer_and_holes<S: Scalar>(
273    outer: &[Vector2<S>],
274    holes: &[Vec<Vector2<S>>],
275) -> Vec<Vector2<S>> {
276    let mut merged: Vec<Vector2<S>> = outer.to_vec();
277    if signed_area2(&merged).definitely_less(S::ZERO) {
278        merged.reverse();
279    }
280
281    for hole in holes {
282        if hole.len() < 3 {
283            continue;
284        }
285        let mut h = hole.clone();
286        if !signed_area2(&h).definitely_less(S::ZERO) {
287            h.reverse();
288        }
289        merge_hole_into(&mut merged, &h);
290    }
291    merged
292}
293
294/// Triangulate `outer` (a simple polygon) minus `holes` (simple polygons
295/// strictly inside `outer` and disjoint from each other), returning
296/// triangles as `(a, b, c)` vertex triples.
297pub fn triangulate_with_holes<S: Scalar>(
298    outer: &[Vector2<S>],
299    holes: &[Vec<Vector2<S>>],
300) -> GeopResult<Vec<(Vector2<S>, Vector2<S>, Vector2<S>)>> {
301    if outer.len() < 3 {
302        return Err(geop_core_math::geop_error::GeopError::new(
303            "triangulate_with_holes: outer polygon needs >= 3 vertices",
304        ));
305    }
306
307    let merged = merge_outer_and_holes(outer, holes);
308    let triples = ear_clip(&merged);
309    Ok(triples
310        .into_iter()
311        .map(|(a, b, c)| (merged[a], merged[b], merged[c]))
312        .collect())
313}
314
315#[cfg(test)]
316mod tests {
317    use super::*;
318    use geop_core_math::for_all_scalars;
319
320    fn check_square_with_hole<S: Scalar>() {
321        let f = S::from_f64;
322        let c = |x: f64, y: f64| Vector2::from_array([f(x), f(y)]);
323        let outer = vec![c(0.0, 0.0), c(1.0, 0.0), c(1.0, 1.0), c(0.0, 1.0)];
324        let hole = vec![c(0.25, 0.25), c(0.75, 0.25), c(0.75, 0.75), c(0.25, 0.75)];
325
326        let tris = triangulate_with_holes(&outer, &[hole]).unwrap();
327        assert!(!tris.is_empty());
328
329        let mut total = S::ZERO;
330        for (a, b, c2) in &tris {
331            total = total.add(cross2(b.sub(a), c2.sub(a)).abs());
332        }
333        // `total` = sum of |2 * signed triangle area| = 2 * (outer area - hole area).
334        let expected = f(2.0 * (1.0 - 0.25));
335        assert!(
336            total.could_be_equal(expected),
337            "total={total}, expected={expected}"
338        );
339    }
340    #[test]
341    fn square_with_hole() {
342        for_all_scalars!(check_square_with_hole);
343    }
344
345    fn check_polygon_signed_area<S: Scalar>() {
346        let f = S::from_f64;
347        let c = |x: f64, y: f64| Vector2::from_array([f(x), f(y)]);
348        let ccw = vec![c(0.0, 0.0), c(1.0, 0.0), c(1.0, 1.0), c(0.0, 1.0)];
349        assert!(polygon_signed_area(&ccw).could_be_equal(f(1.0)));
350
351        let cw: Vec<_> = ccw.into_iter().rev().collect();
352        assert!(polygon_signed_area(&cw).could_be_equal(f(-1.0)));
353    }
354    #[test]
355    fn polygon_signed_area_sign() {
356        for_all_scalars!(check_polygon_signed_area);
357    }
358}
geop_ops_rasterize/polygon_triangulate.rs

geop_ops_rasterize/
polygon_triangulate.rs
