fast_mesh_ops/
bridge.rs

1//! Bridging two boundary loops with a triangle tube.
2//!
3//! [`bridge_loops`] stitches a tube between two boundary loops identified by
4//! [`LoopId`](crate::LoopId). The correspondence is a minimum-area dynamic
5//! program over the loop-by-loop grid with a full start search, so loops of
6//! different lengths zip without twist. The loops must wind in opposite
7//! directions (their Newell vectors must oppose); same-winding loops admit
8//! no outward tube that is manifold at both rims and are rejected. The whole
9//! strip is planned first and committed through one
10//! [`Transaction`](fast_mesh_store::Transaction), so the call is atomic.
11
12use std::cmp::Ordering;
13
14use fast_mesh_store::{Mesh, Transaction, VertexId};
15
16use crate::topology::{BoundaryLoop, EdgeMap, boundary_loops, cross, dot, point, sub};
17use crate::{LoopId, OpsError, Result, push, reserve};
18
19/// Default cap on each loop's length in vertices. Explicit `max_edges`
20/// overrides it; the dynamic program costs O(starts x cells), so raising
21/// the bound past a few hundred gets slow before it gets big.
22pub const BRIDGE_DEFAULT_MAX_EDGES: usize = 128;
23
24/// Options for [`bridge_loops`].
25#[derive(Debug, Clone, Copy, PartialEq)]
26pub struct BridgeOptions {
27    /// The two loops to connect. Order is insignificant: the smaller id is
28    /// always traversed first internally.
29    pub loops: (LoopId, LoopId),
30    /// Maximum vertices per loop. `None` means [`BRIDGE_DEFAULT_MAX_EDGES`].
31    pub max_edges: Option<usize>,
32}
33
34impl BridgeOptions {
35    /// A builder for the two loops to connect.
36    pub fn builder(a: LoopId, b: LoopId) -> BridgeOptionsBuilder {
37        BridgeOptionsBuilder {
38            loops: (a, b),
39            max_edges: None,
40        }
41    }
42}
43
44/// Builder for [`BridgeOptions`].
45#[derive(Debug, Clone, Copy)]
46pub struct BridgeOptionsBuilder {
47    loops: (LoopId, LoopId),
48    max_edges: Option<usize>,
49}
50
51impl BridgeOptionsBuilder {
52    /// Maximum vertices per loop.
53    pub fn max_edges(mut self, max: usize) -> Self {
54        self.max_edges = Some(max);
55        self
56    }
57
58    /// The options as set. Loop validity is checked by [`bridge_loops`].
59    pub fn build(self) -> BridgeOptions {
60        BridgeOptions {
61            loops: self.loops,
62            max_edges: self.max_edges,
63        }
64    }
65}
66
67/// Outcome of [`bridge_loops`].
68#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
69pub struct BridgeReport {
70    /// Triangles added between the loops.
71    pub triangles_added: usize,
72}
73
74/// Stitches a tube between the two loops in `options`.
75///
76/// The first loop (smaller id) is traversed backward and the second
77/// forward; every emitted triangle uses its rim edges opposite to the rim
78/// triangles, so the strip is manifold at both rims, and the minimum-area
79/// correspondence keeps it untwisted. New diagonals must not already exist
80/// in the mesh. Either every triangle lands or nothing does.
81pub fn bridge_loops(tx: &mut Transaction, options: &BridgeOptions) -> Result<BridgeReport> {
82    let (ida, idb) = options.loops;
83    if ida == idb {
84        return Err(OpsError::InvalidInput(
85            "bridge loops must be two distinct loops",
86        ));
87    }
88    let mesh = tx.mesh();
89    let loops = boundary_loops(mesh)?;
90    let la = resolve(&loops, &ida)?;
91    let lb = resolve(&loops, &idb)?;
92    let (first, second) = if loop_key(&ida) <= loop_key(&idb) {
93        (la, lb)
94    } else {
95        (lb, la)
96    };
97    for l in [first, second] {
98        if !l.closed || !l.clean || l.vertices.len() < 3 {
99            return Err(OpsError::InvalidInput(
100                "bridge loops must be closed and clean",
101            ));
102        }
103    }
104    let bound = options.max_edges.unwrap_or(BRIDGE_DEFAULT_MAX_EDGES);
105    if first.vertices.len() > bound || second.vertices.len() > bound {
106        return Err(OpsError::InvalidInput("bridge loop exceeds max_edges"));
107    }
108    for &v in &first.vertices {
109        if second.vertices.contains(&v) {
110            return Err(OpsError::InvalidInput(
111                "bridge loops must not share vertices",
112            ));
113        }
114    }
115    for l in [first, second] {
116        for &v in &l.vertices {
117            if !point(mesh, v).iter().all(|c| c.is_finite()) {
118                return Err(OpsError::InvalidInput(
119                    "bridge loop positions must be finite",
120                ));
121            }
122        }
123    }
124    let na = newell(mesh, &first.vertices);
125    let nb = newell(mesh, &second.vertices);
126    if !matches!(dot(na, nb).partial_cmp(&0.0), Some(Ordering::Less)) {
127        return Err(OpsError::InvalidInput(
128            "bridge loops must wind in opposite directions",
129        ));
130    }
131    let strip = best_strip(mesh, &first.vertices, &second.vertices)?;
132    for t in &strip {
133        let area2 = tri_area2(&point(mesh, t[0]), &point(mesh, t[1]), &point(mesh, t[2]));
134        if area2 == 0.0 {
135            return Err(OpsError::InvalidInput("bridge triangulation is degenerate"));
136        }
137    }
138    let edges = EdgeMap::build(mesh)?;
139    let mut chords = chords_of(&strip, &first.vertices)?;
140    chords.sort_unstable();
141    chords.dedup();
142    for &(a, b) in &chords {
143        if edges.kind(a, b).is_some() {
144            return Err(OpsError::InvalidInput("bridge diagonal already exists"));
145        }
146    }
147    let report = BridgeReport {
148        triangles_added: strip.len(),
149    };
150    tx.add_triangles(&strip)?;
151    Ok(report)
152}
153
154fn loop_key(id: &LoopId) -> (u32, u32, u64) {
155    (id.min_vertex, id.len, id.hash)
156}
157
158/// Finds a loop by id in a [`boundary_loops`] vec, like fill does.
159fn resolve<'a>(loops: &'a [BoundaryLoop], id: &LoopId) -> Result<&'a BoundaryLoop> {
160    let at = loops.partition_point(|l| l.id.min_vertex < id.min_vertex);
161    loops
162        .get(at)
163        .filter(|l| l.id == *id)
164        .ok_or(OpsError::StaleLoops)
165}
166
167/// The Newell vector of a loop: its direction is the loop's winding normal.
168fn newell(mesh: &Mesh, verts: &[VertexId]) -> [f64; 3] {
169    let mut n = [0.0; 3];
170    for k in 0..verts.len() {
171        let c = cross(
172            point(mesh, verts[k]),
173            point(mesh, verts[(k + 1) % verts.len()]),
174        );
175        n[0] += c[0];
176        n[1] += c[1];
177        n[2] += c[2];
178    }
179    n
180}
181
182/// Twice the area of a triangle over f64 positions.
183fn tri_area2(a: &[f64; 3], b: &[f64; 3], c: &[f64; 3]) -> f64 {
184    let n = cross(sub(*b, *a), sub(*c, *a));
185    dot(n, n).sqrt()
186}
187
188/// The minimum-area strip over all start offsets of the second loop. The
189/// first loop is traversed backward, the second forward; see the
190/// [`bridge_loops`] docs for why that combination is manifold.
191fn best_strip(mesh: &Mesh, first: &[VertexId], second: &[VertexId]) -> Result<Vec<[VertexId; 3]>> {
192    let n = first.len();
193    let m = second.len();
194    // X walks the first loop backward from its start vertex; Y walks the
195    // second loop forward, doubled so every start offset is a plain slice.
196    let mut xa: Vec<VertexId> = Vec::new();
197    reserve(&mut xa, n + 1)?;
198    for i in 0..=n {
199        push(&mut xa, first[(n - i) % n])?;
200    }
201    let mut px: Vec<[f64; 3]> = Vec::new();
202    reserve(&mut px, n + 1)?;
203    for &v in &xa {
204        push(&mut px, point(mesh, v))?;
205    }
206    let mut py: Vec<[f64; 3]> = Vec::new();
207    reserve(&mut py, 2 * m + 1)?;
208    for k in 0..=2 * m {
209        push(&mut py, point(mesh, second[k % m]))?;
210    }
211    let cells = (n + 1) * (m + 1);
212    let mut cost: Vec<f64> = Vec::new();
213    reserve(&mut cost, cells)?;
214    cost.extend(std::iter::repeat_n(f64::INFINITY, cells));
215    let mut back: Vec<u8> = Vec::new();
216    reserve(&mut back, cells)?;
217    back.extend(std::iter::repeat_n(0, cells));
218    let mut best = f64::INFINITY;
219    let mut best_s = 0;
220    for s in 0..m {
221        let total = dp_fill(&px, &py[s..s + m + 1], &mut cost, &mut back);
222        if total < best {
223            best = total;
224            best_s = s;
225        }
226    }
227    dp_fill(&px, &py[best_s..best_s + m + 1], &mut cost, &mut back);
228    // Backtrace the staircase. Rim states are forced along their only
229    // incoming edge so garbage table entries cannot underflow the walk.
230    let mut strip: Vec<[VertexId; 3]> = Vec::new();
231    reserve(&mut strip, n + m)?;
232    let (mut i, mut k) = (n, m);
233    while i > 0 || k > 0 {
234        let step_a = if k == 0 {
235            true
236        } else if i == 0 {
237            false
238        } else {
239            back[i * (m + 1) + k] == 0
240        };
241        if step_a {
242            push(&mut strip, [xa[i - 1], xa[i], second[(best_s + k) % m]])?;
243            i -= 1;
244        } else {
245            push(
246                &mut strip,
247                [
248                    xa[i],
249                    second[(best_s + k) % m],
250                    second[(best_s + k - 1) % m],
251                ],
252            )?;
253            k -= 1;
254        }
255    }
256    strip.reverse();
257    Ok(strip)
258}
259
260/// Fills the correspondence table for one start offset and returns the
261/// minimum total doubled area. Moving along X emits (x, x', y); moving
262/// along Y emits (x, y', y); ties prefer X. Areas only, no ids here.
263fn dp_fill(px: &[[f64; 3]], yb: &[[f64; 3]], cost: &mut [f64], back: &mut [u8]) -> f64 {
264    let n = px.len() - 1;
265    let m = yb.len() - 1;
266    for c in cost.iter_mut() {
267        *c = f64::INFINITY;
268    }
269    cost[0] = 0.0;
270    for i in 0..=n {
271        for k in 0..=m {
272            if i == 0 && k == 0 {
273                continue;
274            }
275            let mut best = f64::INFINITY;
276            let mut dir = 0u8;
277            if i > 0 {
278                let c = cost[(i - 1) * (m + 1) + k] + tri_area2(&px[i - 1], &px[i], &yb[k]);
279                if c <= best {
280                    best = c;
281                    dir = 0;
282                }
283            }
284            if k > 0 {
285                let c = cost[i * (m + 1) + k - 1] + tri_area2(&px[i], &yb[k], &yb[k - 1]);
286                if c < best {
287                    best = c;
288                    dir = 1;
289                }
290            }
291            cost[i * (m + 1) + k] = best;
292            back[i * (m + 1) + k] = dir;
293        }
294    }
295    cost[n * (m + 1) + m]
296}
297
298/// The strip's mixed first/second-loop edges: rim edges stay within one
299/// loop, so any edge the strip shares with the mesh must be one of these.
300/// Every strip vertex belongs to a loop, so "not first" means "second".
301fn chords_of(strip: &[[VertexId; 3]], first: &[VertexId]) -> Result<Vec<(VertexId, VertexId)>> {
302    let mut a_sorted: Vec<VertexId> = Vec::new();
303    reserve(&mut a_sorted, first.len())?;
304    a_sorted.extend_from_slice(first);
305    a_sorted.sort_unstable();
306    let mut chords: Vec<(VertexId, VertexId)> = Vec::new();
307    reserve(&mut chords, 2 * strip.len())?;
308    for t in strip {
309        for (p, q) in [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])] {
310            let pa = a_sorted.binary_search(&p).is_ok();
311            let qa = a_sorted.binary_search(&q).is_ok();
312            if pa != qa {
313                push(&mut chords, (p.min(q), p.max(q)))?;
314            }
315        }
316    }
317    Ok(chords)
318}
319
320#[cfg(test)]
321mod tests {
322    use fast_mesh_store::{Mesh, Transaction};
323
324    use super::*;
325    use crate::test_util::{edit, mesh_of, parts, sweep};
326    use crate::topology::EdgeMap;
327    use crate::{OpsError, analyze, boundary_loops};
328
329    /// Quad A at z=0 facing -z, quad B at z=1 facing +z: both rims face
330    /// away from the gap, so an outward tube is manifold at both rims.
331    fn two_quads() -> Mesh {
332        mesh_of(
333            &[
334                0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, //
335                0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0,
336            ],
337            &[[0, 2, 1], [0, 3, 2], [4, 5, 6], [4, 6, 7]],
338        )
339    }
340
341    /// Triangle at z=0 facing -z, quad at z=1 facing +z.
342    fn tri_square() -> Mesh {
343        mesh_of(
344            &[
345                0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, //
346                0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0,
347            ],
348            &[[0, 2, 1], [3, 4, 5], [3, 5, 6]],
349        )
350    }
351
352    fn loop_of(mesh: &Mesh, min_vertex: u32) -> crate::LoopId {
353        boundary_loops(mesh)
354            .unwrap()
355            .iter()
356            .find(|l| l.id.min_vertex == min_vertex)
357            .unwrap()
358            .id
359    }
360
361    /// Every triangle added past `before` faces away from the tube axis.
362    fn assert_outward(mesh: &Mesh, before: usize) {
363        let (pos, tris) = parts(mesh);
364        let p = |v: u32| -> [f64; 3] {
365            let i = 3 * v as usize;
366            [
367                f64::from(pos[i]),
368                f64::from(pos[i + 1]),
369                f64::from(pos[i + 2]),
370            ]
371        };
372        // Axis through the loop centroids: A is verts 0..4, B the rest.
373        let mut ca = [0.0; 3];
374        let mut cb = [0.0; 3];
375        for v in 0..4u32 {
376            let q = p(v);
377            for k in 0..3 {
378                ca[k] += q[k] / 4.0;
379            }
380        }
381        for v in 4..8u32 {
382            let q = p(v);
383            for k in 0..3 {
384                cb[k] += q[k] / 4.0;
385            }
386        }
387        assert!(tris.len() > before);
388        for t in &tris[before..] {
389            let (a, b, c) = (p(t[0]), p(t[1]), p(t[2]));
390            let n = [
391                (b[1] - a[1]) * (c[2] - a[2]) - (b[2] - a[2]) * (c[1] - a[1]),
392                (b[2] - a[2]) * (c[0] - a[0]) - (b[0] - a[0]) * (c[2] - a[2]),
393                (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]),
394            ];
395            let m = [
396                (a[0] + b[0] + c[0]) / 3.0 - (ca[0] + cb[0]) / 2.0,
397                (a[1] + b[1] + c[1]) / 3.0 - (ca[1] + cb[1]) / 2.0,
398                0.0,
399            ];
400            let d = n[0] * m[0] + n[1] * m[1] + n[2] * m[2];
401            assert!(d > 0.0, "inward or degenerate: {t:?} d={d}");
402        }
403    }
404
405    #[test]
406    fn bridge_zips_two_square_loops() {
407        let mut mesh = two_quads();
408        let before = parts(&mesh).1.len();
409        let a = loop_of(&mesh, 0);
410        let b = loop_of(&mesh, 4);
411        let options = BridgeOptions::builder(a, b).build();
412        let (report, _) = edit(&mut mesh, |tx| bridge_loops(tx, &options));
413        assert_eq!(report.triangles_added, 8);
414        assert_outward(&mesh, before);
415        // Every rim edge is interior now, and so is every new diagonal.
416        let edges = EdgeMap::build(&mesh).unwrap();
417        let (_, tris) = parts(&mesh);
418        for t in &tris[before..] {
419            for (x, y) in [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])] {
420                assert_eq!(edges.triangles(x, y).len(), 2, "edge ({x}, {y})");
421            }
422        }
423        let conflicts = analyze(&mesh, &Default::default())
424            .unwrap()
425            .orientation_conflicts;
426        assert_eq!(conflicts, 0);
427    }
428
429    #[test]
430    fn bridge_handles_mismatched_loop_sizes() {
431        let mut mesh = tri_square();
432        let before = parts(&mesh).1.len();
433        let a = loop_of(&mesh, 0);
434        let b = loop_of(&mesh, 3);
435        let options = BridgeOptions::builder(a, b).build();
436        let (report, _) = edit(&mut mesh, |tx| bridge_loops(tx, &options));
437        assert_eq!(report.triangles_added, 7);
438        let edges = EdgeMap::build(&mesh).unwrap();
439        let (_, tris) = parts(&mesh);
440        for t in &tris[before..] {
441            for (x, y) in [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])] {
442                assert_eq!(edges.triangles(x, y).len(), 2, "edge ({x}, {y})");
443            }
444        }
445        let conflicts = analyze(&mesh, &Default::default())
446            .unwrap()
447            .orientation_conflicts;
448        assert_eq!(conflicts, 0);
449    }
450
451    #[test]
452    fn bridge_rejects_same_winding() {
453        // Both quads face +z (both rims face into the gap on one side):
454        // no outward tube can be manifold at both rims.
455        let mut mesh = mesh_of(
456            &[
457                0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, //
458                0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0,
459            ],
460            &[[0, 1, 2], [0, 2, 3], [4, 5, 6], [4, 6, 7]],
461        );
462        let before = parts(&mesh);
463        let a = loop_of(&mesh, 0);
464        let b = loop_of(&mesh, 4);
465        let options = BridgeOptions::builder(a, b).build();
466        let mut tx = mesh.begin_edit();
467        assert!(matches!(
468            bridge_loops(&mut tx, &options),
469            Err(OpsError::InvalidInput(_))
470        ));
471        drop(tx);
472        assert_eq!(parts(&mesh), before);
473    }
474
475    #[test]
476    fn bridge_rejects_bad_selection() {
477        let mut mesh = two_quads();
478        let before = parts(&mesh);
479        let a = loop_of(&mesh, 0);
480        let b = loop_of(&mesh, 4);
481        // Same loop twice.
482        let options = BridgeOptions::builder(a, a).build();
483        let mut tx = mesh.begin_edit();
484        assert!(matches!(
485            bridge_loops(&mut tx, &options),
486            Err(OpsError::InvalidInput(_))
487        ));
488        drop(tx);
489        // Unknown loop id.
490        let fake = crate::LoopId {
491            min_vertex: 99,
492            len: 4,
493            hash: 7,
494        };
495        let options = BridgeOptions::builder(a, fake).build();
496        let mut tx = mesh.begin_edit();
497        assert_eq!(bridge_loops(&mut tx, &options), Err(OpsError::StaleLoops));
498        drop(tx);
499        // Over the per-loop bound.
500        let options = BridgeOptions::builder(a, b).max_edges(3).build();
501        let mut tx = mesh.begin_edit();
502        assert!(matches!(
503            bridge_loops(&mut tx, &options),
504            Err(OpsError::InvalidInput(_))
505        ));
506        drop(tx);
507        assert_eq!(parts(&mesh), before);
508    }
509
510    #[test]
511    fn bridge_sweep_keeps_atomicity() {
512        let mesh = two_quads();
513        let a = loop_of(&mesh, 0);
514        let b = loop_of(&mesh, 4);
515        let options = BridgeOptions::builder(a, b).build();
516        let n = sweep(mesh, |tx: &mut Transaction| {
517            bridge_loops(tx, &options)?;
518            Ok(())
519        });
520        assert!(n > 0);
521    }
522}
fast_mesh_ops/bridge.rs

fast_mesh_ops/
bridge.rs
