axiolid_triangulate/
build.rs

1// SPDX-License-Identifier: MPL-2.0
2
3//! Incremental constrained Delaunay construction.
4
5use axiolid_core::Point2;
6use axiolid_guarantees::Sign;
7use axiolid_predicates::orient2d;
8
9use crate::mesh::{Triangulation, TriangulationError, NO_HALFEDGE};
10use crate::recover::recover_constraints;
11use crate::{decided, in_circumcircle, turns_left, Constraint};
12
13// ---------------------------------------------------------------------------
14// Construction
15// ---------------------------------------------------------------------------
16
17/// Build a constrained Delaunay triangulation of `points`.
18///
19/// Every edge in `constraints` appears in the output as an edge of some
20/// triangle. Away from constraints the result satisfies the Delaunay
21/// empty-circumcircle property, decided exactly.
22///
23/// # Errors
24///
25/// Returns [`TriangulationError`] when the input has fewer than three points,
26/// is entirely collinear, references a vertex that does not exist, or
27/// contains constraints that cross one another.
28pub fn triangulate(
29    points: &[Point2],
30    constraints: &[Constraint],
31) -> Result<Triangulation, TriangulationError> {
32    if points.len() < 3 {
33        return Err(TriangulationError::TooFewPoints);
34    }
35    let count = u32::try_from(points.len()).unwrap_or(u32::MAX);
36    for c in constraints {
37        if c.a >= count || c.b >= count {
38            let index = if c.a >= count { c.a } else { c.b };
39            return Err(TriangulationError::ConstraintOutOfRange { index });
40        }
41    }
42    if !has_non_collinear_triple(points) {
43        return Err(TriangulationError::AllPointsCollinear);
44    }
45
46    let mut state = Builder::new(points);
47    state.insert_all()?;
48    let mut triangulation = state.finish(points, constraints);
49    recover_constraints(&mut triangulation)?;
50    Ok(triangulation)
51}
52
53/// Whether the input contains three points that are not collinear.
54///
55/// Checked with the exact predicate rather than by an area threshold: a
56/// tolerance here would reject a legitimately thin but non-degenerate input.
57fn has_non_collinear_triple(points: &[Point2]) -> bool {
58    let a = points[0];
59    // Find the first point distinct from `a`, then the first that is not on
60    // the line through the two. Scanning is linear and runs once.
61    let Some(b) = points.iter().copied().find(|p| *p != a) else {
62        return false;
63    };
64    points
65        .iter()
66        .any(|&c| turns_left(a, b, c) || turns_left(b, a, c))
67}
68
69/// Incremental Delaunay construction over a super triangle.
70struct Builder {
71    points: Vec<Point2>,
72    triangles: Vec<u32>,
73    halfedges: Vec<u32>,
74    real_count: usize,
75}
76
77impl Builder {
78    fn new(points: &[Point2]) -> Self {
79        let mut all = points.to_vec();
80        // A super triangle that strictly contains every input point. Scaled
81        // generously off the bounding box: a tight super triangle puts input
82        // points on its edges, and a point exactly on an edge has no
83        // containing triangle to split.
84        let (min, max) = bounds(points);
85        let dx = max.x - min.x;
86        let dy = max.y - min.y;
87        let span = if dx > dy { dx } else { dy };
88        let span = if span > 0.0 { span } else { 1.0 };
89        let cx = (min.x + max.x) * 0.5;
90        let cy = (min.y + max.y) * 0.5;
91        let far = span * 32.0;
92        all.push(Point2::new(cx - far, cy - far));
93        all.push(Point2::new(cx + far, cy - far));
94        all.push(Point2::new(cx, cy + far));
95
96        let n = points.len() as u32;
97        Self {
98            points: all,
99            triangles: vec![n, n + 1, n + 2],
100            halfedges: vec![NO_HALFEDGE, NO_HALFEDGE, NO_HALFEDGE],
101            real_count: points.len(),
102        }
103    }
104
105    fn insert_all(&mut self) -> Result<(), TriangulationError> {
106        for v in 0..self.real_count as u32 {
107            self.insert_point(v)?;
108        }
109        Ok(())
110    }
111
112    /// Insert one vertex, splitting its containing triangle and restoring the
113    /// Delaunay property by flipping.
114    fn insert_point(&mut self, v: u32) -> Result<(), TriangulationError> {
115        // A vertex that cannot be located is a corrupted adjacency, not a
116        // benign skip. Dropping it silently produced a triangulation missing
117        // an input point -- the caller would get a valid-looking mesh with a
118        // hole where their vertex should be.
119        let Some(t) = self.locate(self.points[v as usize]) else {
120            return Err(TriangulationError::VertexUnplaceable { index: v });
121        };
122        let (a, b, c) = (
123            self.triangles[3 * t],
124            self.triangles[3 * t + 1],
125            self.triangles[3 * t + 2],
126        );
127        let p = self.points[v as usize];
128
129        // On an edge: split the edge, and the triangles either side of it,
130        // rather than the one triangle into three -- one of the three would
131        // have no area, and flips over a triangle with no area never settle.
132        for i in 0..3 {
133            let (u, w) = (
134                self.triangles[3 * t + i],
135                self.triangles[3 * t + (i + 1) % 3],
136            );
137            if decided(orient2d(
138                self.points[u as usize],
139                self.points[w as usize],
140                p,
141            )) != Sign::Zero
142            {
143                continue;
144            }
145            let twin = self.halfedges[3 * t + i];
146            if twin == NO_HALFEDGE {
147                return Err(TriangulationError::VertexUnplaceable { index: v });
148            }
149            let x = self.triangles[3 * t + (i + 2) % 3];
150            let ot = twin as usize / 3;
151            let oi = twin as usize % 3;
152            let y = self.triangles[3 * ot + (oi + 2) % 3];
153            // (u, w, x) and (w, u, y) become (u, v, x), (v, w, x),
154            // (w, v, y), (v, u, y).
155            let t1 = self.triangles.len() / 3;
156            let t2 = t1 + 1;
157            self.triangles[3 * t] = x;
158            self.triangles[3 * t + 1] = u;
159            self.triangles[3 * t + 2] = v;
160            self.triangles[3 * ot] = y;
161            self.triangles[3 * ot + 1] = w;
162            self.triangles[3 * ot + 2] = v;
163            self.triangles.extend_from_slice(&[w, x, v]);
164            self.triangles.extend_from_slice(&[u, y, v]);
165            self.halfedges.extend_from_slice(&[NO_HALFEDGE; 6]);
166            self.rebuild_adjacency();
167            for tri in [t, ot, t1, t2] {
168                self.legalize(3 * tri);
169            }
170            return Ok(());
171        }
172
173        // Reuse slot `t` for the first sub-triangle and append the other two,
174        // so existing neighbour indices into `t` stay valid where possible.
175        let t1 = self.triangles.len() / 3;
176        let t2 = t1 + 1;
177        self.triangles[3 * t] = a;
178        self.triangles[3 * t + 1] = b;
179        self.triangles[3 * t + 2] = v;
180        self.triangles.extend_from_slice(&[b, c, v]);
181        self.triangles.extend_from_slice(&[c, a, v]);
182        self.halfedges.extend_from_slice(&[NO_HALFEDGE; 6]);
183
184        // Rebuild adjacency from the triangle list rather than patching the
185        // six affected twin pointers by hand. Hand-patching was wrong here:
186        // reusing slot `t` for one sub-triangle leaves the OLD neighbours of
187        // `t` pointing at halfedges that now belong to a different triangle,
188        // and the corruption only surfaces later as a walk that leaves the
189        // mesh. The rebuild is O(n) per insertion, which is the price of a
190        // provably consistent structure; `locate` stays correct, and the
191        // refinement loop above is what dominates runtime anyway.
192        self.rebuild_adjacency();
193
194        self.legalize(3 * t);
195        self.legalize(3 * t1);
196        self.legalize(3 * t2);
197        Ok(())
198    }
199
200    /// Recompute every twin pointer from the triangle list.
201    fn rebuild_adjacency(&mut self) {
202        use std::collections::HashMap;
203        let count = self.triangles.len();
204        self.halfedges.clear();
205        self.halfedges.resize(count, NO_HALFEDGE);
206        let mut seen: HashMap<(u32, u32), u32> = HashMap::with_capacity(count);
207        for e in 0..count {
208            let t = e / 3;
209            let i = e % 3;
210            let from = self.triangles[3 * t + i];
211            let to = self.triangles[3 * t + (i + 1) % 3];
212            if let Some(&twin) = seen.get(&(to, from)) {
213                self.halfedges[e] = twin;
214                self.halfedges[twin as usize] = e as u32;
215            } else {
216                seen.insert((from, to), e as u32);
217            }
218        }
219    }
220
221    /// Find a triangle containing `p` by walking from the last triangle.
222    ///
223    /// A straight walk is O(sqrt n) on well-distributed input and needs no
224    /// auxiliary structure. Returns `None` only if the walk leaves the mesh,
225    /// which the super triangle makes impossible for in-range points.
226    fn locate(&self, p: Point2) -> Option<usize> {
227        let mut t = self.triangles.len() / 3 - 1;
228        for _ in 0..self.triangles.len() {
229            let (a, b, c) = (
230                self.points[self.triangles[3 * t] as usize],
231                self.points[self.triangles[3 * t + 1] as usize],
232                self.points[self.triangles[3 * t + 2] as usize],
233            );
234            // Step across the first edge that `p` lies strictly outside of.
235            // For a counter-clockwise triangle, "outside edge (u, w)" means
236            // p is strictly to the RIGHT of u->w. Testing `turns_left(w, u, p)`
237            // instead also fires when p is exactly ON the edge, which sends
238            // the walk across a boundary it should have stopped at.
239            let mut moved = false;
240            for (i, (u, w)) in [(a, b), (b, c), (c, a)].into_iter().enumerate() {
241                let right_of = decided(orient2d(u, w, p)) == Sign::Negative;
242                if right_of {
243                    let twin = self.halfedges[3 * t + i];
244                    if twin == NO_HALFEDGE {
245                        return None;
246                    }
247                    t = twin as usize / 3;
248                    moved = true;
249                    break;
250                }
251            }
252            if !moved {
253                return Some(t);
254            }
255        }
256        None
257    }
258
259    /// Restore the Delaunay property across `edge`, whose triangle's third
260    /// corner is the vertex just inserted, and propagate to the edges
261    /// across from it.
262    fn legalize(&mut self, edge: usize) {
263        let mut stack = vec![edge];
264        // Lawson's flips after one insertion are few; the bound only keeps
265        // a corrupted mesh from spinning. Anything left is repaired after
266        // constraint recovery.
267        let mut budget = 4 * self.triangles.len() + 64;
268        while let Some(e) = stack.pop() {
269            if budget == 0 {
270                break;
271            }
272            budget -= 1;
273            let twin = self.halfedges[e];
274            if twin == NO_HALFEDGE {
275                continue;
276            }
277            let t = e / 3;
278            let ti = e % 3;
279            let o = twin as usize;
280            let ot = o / 3;
281            let oi = o % 3;
282
283            let p0 = self.triangles[3 * t + ti];
284            let p1 = self.triangles[3 * t + (ti + 1) % 3];
285            let apex = self.triangles[3 * t + (ti + 2) % 3];
286            let other = self.triangles[3 * ot + (oi + 2) % 3];
287
288            if !in_circumcircle(
289                self.points[p0 as usize],
290                self.points[p1 as usize],
291                self.points[apex as usize],
292                self.points[other as usize],
293            ) {
294                continue;
295            }
296            // Only a strictly convex quadrilateral flips without inverting.
297            let q = |k: u32| self.points[k as usize];
298            if !(turns_left(q(apex), q(p0), q(other)) && turns_left(q(other), q(p1), q(apex))) {
299                continue;
300            }
301
302            // Flip the shared edge to the other diagonal: (p0, other, apex)
303            // and (p1, apex, other).
304            self.triangles[3 * t + (ti + 1) % 3] = other;
305            self.triangles[3 * ot + (oi + 1) % 3] = apex;
306
307            // Same reasoning as insertion: recompute adjacency rather than
308            // re-point twins by hand.
309            self.rebuild_adjacency();
310
311            // The two edges now across from the inserted vertex: p0 -> other
312            // and other -> p1. (Pushing the new diagonal instead, as this
313            // did, checked an edge that is Delaunay by construction and left
314            // these unchecked; thin quads then flipped back and forth.)
315            stack.push(3 * t + ti);
316            stack.push(3 * ot + (oi + 2) % 3);
317        }
318    }
319
320    /// Drop the super triangle and everything touching it.
321    fn finish(self, original: &[Point2], constraints: &[Constraint]) -> Triangulation {
322        let limit = self.real_count as u32;
323        let mut triangles = Vec::with_capacity(self.triangles.len());
324        for t in self.triangles.chunks_exact(3) {
325            if t.iter().all(|&v| v < limit) {
326                triangles.extend_from_slice(t);
327            }
328        }
329        let mut sorted = constraints.to_vec();
330        sorted.sort_unstable();
331        sorted.dedup();
332        let mut out = Triangulation {
333            points: original.to_vec(),
334            triangles,
335            halfedges: Vec::new(),
336            constraints: sorted,
337        };
338        out.rebuild_halfedges();
339        out
340    }
341}
342
343/// Axis-aligned bounds of a point set.
344fn bounds(points: &[Point2]) -> (Point2, Point2) {
345    let mut min = points[0];
346    let mut max = points[0];
347    for p in &points[1..] {
348        min.x = min.x.min(p.x);
349        min.y = min.y.min(p.y);
350        max.x = max.x.max(p.x);
351        max.y = max.y.max(p.y);
352    }
353    (min, max)
354}
axiolid_triangulate/build.rs

axiolid_triangulate/
build.rs
