From e3459edea80685c26d53c8b60165e8e786d224c3 Mon Sep 17 00:00:00 2001 From: Irakli Kourtis Date: Wed, 14 Oct 2026 15:01:05 -0400 Subject: [PATCH] Display tessellation: trimmed faces as four-sided patches (#303) (#1163) Trimmed faces are cut into four-sided patches (squares around holes, butterflies for round faces, motorcycle traces from concave corners) and each is meshed as a Coons grid, with a cut grid as the fallback; cells are at most 2:1 along straight directions. Folded blends are refused rather than overlapped. The README lists display meshes and the gpu crate. --- README.md | 12 +- crates/kernel/src/display_mesh_tests.rs | 4 +- crates/kernel/src/tess/cdt.rs | 63 +++ crates/kernel/src/tess/collar.rs | 224 ++++++++++ crates/kernel/src/tess/coons.rs | 306 +++++++++++++- crates/kernel/src/tess/cutgrid.rs | 518 ++++++++++++++++++++++++ crates/kernel/src/tess/mod.rs | 240 ++++++++++- crates/kernel/src/tess/motorcycle.rs | 354 ++++++++++++++++ crates/kernel/src/tess/pieces.rs | 348 ++++++++++++++++ crates/kernel/src/tess/size.rs | 16 +- crates/kernel/src/tess/tests.rs | 21 +- 11 files changed, 2068 insertions(+), 38 deletions(-) create mode 100644 crates/kernel/src/tess/collar.rs create mode 100644 crates/kernel/src/tess/cutgrid.rs create mode 100644 crates/kernel/src/tess/motorcycle.rs create mode 100644 crates/kernel/src/tess/pieces.rs diff --git a/README.md b/README.md index 815cc4d3..c9990f2c 100644 --- a/README.md +++ b/README.md @@ -134,6 +134,16 @@ rebuild with no errors, replayed here at 2× speed. (manifold solid B-reps, or trimmed and bounded surfaces sewn into solids) as a base feature that later features build on — open a `.step`/`.stp`/`.igs`/`.iges` file, insert one into a design, or drop it on the window ([docs/step-import.md](docs/step-import.md)). +- **Display meshes** from our own tessellation planner: each face becomes a plan of (u, v) + points and triangles that shares its edge nodes with its neighbours, so the mesh is watertight. + Facets follow the chord tolerance, the facet angle and a longest side, with cells at most 2:1 + along straight directions. Box-bounded and four-cornered faces mesh as structured grids. + Trimmed faces are cut into four-sided patches: squares around holes, butterflies in round + faces, motorcycle traces from concave corners (after Lyon et al. 2019). Each patch then meshes + as a Coons grid. Intersection edges are evaluated by Newton on both surfaces. A compute shader + evaluates the plans for the whole scene in one dispatch, with a CPU fallback when the GPU + can't build it. Planning and evaluating the 78 reference bodies takes about 65 ms, 50–65× + faster than the previous mesher. - **Desktop app** (egui + wgpu): toolbar with workspace tabs, browser, timeline, 3D viewport with axis view and navigation bar, sketch tools, feature dialogs, command palette, and a JSON-lines control channel so agents and tests can drive it. @@ -215,7 +225,7 @@ configuration and an example session: [docs/mcp.md](docs/mcp.md). | L2 | `doc` | parameters, expressions, the feature timeline and its incremental evaluation | | L3 | `io` | design files, STL/OBJ/STEP/3MF export, STEP and 3MF/STL import features | | L4 | `engine` | the session and the command registry (everything is a command) | -| L5 | `ui-egui`, `mcp` | the swappable desktop front end; the MCP server (headless or bridged to the app) | +| L5 | `gpu`, `ui-egui`, `mcp` | compute-shader evaluation of display tessellation plans; the swappable desktop front end; the MCP server (headless or bridged to the app) | | apps | `solidcraft`, `solidcraft-cli` | desktop app, headless CLI | diff --git a/crates/kernel/src/display_mesh_tests.rs b/crates/kernel/src/display_mesh_tests.rs index f2b0a17f..13203e9f 100644 --- a/crates/kernel/src/display_mesh_tests.rs +++ b/crates/kernel/src/display_mesh_tests.rs @@ -337,7 +337,9 @@ fn round_faces_of_every_kind_are_divided_by_the_angle() { let check = |name: &str, b: &Body| { let coarse = b.display_mesh(3.0).unwrap(); let round = b.display_mesh_angular(3.0, DISPLAY_ANGLE).unwrap(); - assert!(turn(b, &coarse) > 2.0 * DISPLAY_ANGLE, "{name}: the chord tolerance leaves facets of {} degrees", turn(b, &coarse).to_degrees()); + // The chord tolerance alone leaves facets coarser than the angle (the size rules' longest + // side and aspect limit already keep them near 40°), so the angle has work to do. + assert!(turn(b, &coarse) > 1.5 * DISPLAY_ANGLE, "{name}: the chord tolerance leaves facets of {} degrees", turn(b, &coarse).to_degrees()); // A side crosses at most one facet in each direction. assert!(turn(b, &round) <= 1.5 * DISPLAY_ANGLE, "{name}: a side turns {} degrees", turn(b, &round).to_degrees()); assert_eq!(open_sides(&round, 1e-6), 0, "{name}"); diff --git a/crates/kernel/src/tess/cdt.rs b/crates/kernel/src/tess/cdt.rs index 1dfed9c0..ec0df6d4 100644 --- a/crates/kernel/src/tess/cdt.rs +++ b/crates/kernel/src/tess/cdt.rs @@ -135,6 +135,69 @@ pub(crate) fn triangulate(boundary: &[[f64; 2]], loops: &[std::ops::Range Some((extra, tris)) } +/// The band between a face's boundary `loops` (ranges of the first `n` points of `points`) and +/// the staircase around a block of grid cells inside it (`stairs`, sides between grid points +/// of `points`): a constrained Delaunay triangulation of both with no point added, keeping the +/// faces inside the boundary and outside the block (an odd number of constraints from the +/// outside). Distances in the metric `scale`. Triangles counter-clockwise in (u, v), indices +/// into `points`; `None` when a constraint cannot be kept (the staircase meets the boundary). +pub(crate) fn band(points: &[[f64; 2]], n: usize, loops: &[std::ops::Range], stairs: &[(u32, u32)], scale: [f64; 2]) -> Option> { + let to = |p: &[f64; 2]| spade::mitigate_underflow(Point2::new(p[0] * scale[0], p[1] * scale[1])); + let mut cdt = Cdt::new(); + let mut index: HashMap = HashMap::new(); + let mut handle: HashMap = HashMap::new(); + let mut add = |cdt: &mut Cdt, i: u32| -> Option { + if let Some(h) = handle.get(&i) { + return Some(*h); + } + let p = points.get(i as usize)?; + if !p.iter().all(|x| x.is_finite()) { + return None; + } + let h = cdt.insert(to(p)).ok()?; + // Two points at one place (pole copies) share a vertex: the first index names it. + index.entry(h).or_insert(i); + handle.insert(i, h); + Some(h) + }; + let mut sides: Vec<(FixedVertexHandle, FixedVertexHandle)> = Vec::new(); + for r in loops { + let len = r.len(); + for k in 0..len { + let (a, b) = (u32::try_from(r.start + k).ok()?, u32::try_from(r.start + (k + 1) % len).ok()?); + if (a as usize) >= n || (b as usize) >= n { + return None; + } + sides.push((add(&mut cdt, a)?, add(&mut cdt, b)?)); + } + } + for &(a, b) in stairs { + sides.push((add(&mut cdt, a)?, add(&mut cdt, b)?)); + } + for (a, b) in sides { + if a == b || cdt.exists_constraint(a, b) { + continue; + } + if !cdt.can_add_constraint(a, b) { + return None; + } + cdt.add_constraint(a, b); + } + let region = region(&cdt); + let mut tris = Vec::new(); + for face in cdt.inner_faces() { + if !region.get(face.index()).copied().unwrap_or(false) { + continue; + } + let mut t = [0u32; 3]; + for (slot, v) in t.iter_mut().zip(face.vertices()) { + *slot = *index.get(&v.fix())?; + } + tris.push(t); + } + Some(tris) +} + /// Insert the parameters of the surface point under every triangle centre that lies more than /// the tolerance from the triangle or whose corner normals turn more than the facet angle, and /// under the middle of every non-constraint side that sags more than the tolerance or is diff --git a/crates/kernel/src/tess/collar.rs b/crates/kernel/src/tess/collar.rs new file mode 100644 index 00000000..341271b8 --- /dev/null +++ b/crates/kernel/src/tess/collar.rs @@ -0,0 +1,224 @@ +//! The cuts that make every patch of a trimmed face four-sided. Each hole gets a square around +//! it in (u, v) (sides on iso-lines, halfway to the nearest other feature), joined to the hole at +//! the hole's four 45° points: the ring between hole and square is four collar patches (a quarter +//! of the hole's edge, a straight cut, a side of the square, a straight cut), the O-grid round a +//! hole. The motorcycle graph then runs with the squares as the obstacles (Lyon et al. 2019 §4): +//! a square's corners are its concave peaks, so the traces extend the squares' sides, and every +//! patch outside the collars is a rectangle bounded by iso-lines. + +use super::motorcycle::{self, End as Stop}; +use super::pieces::{End, Segment}; + +/// The face's cuts, given its boundary loops (ranges of `uvs`), the metric `scale` and the face's +/// seed (its node spacing, mm; traces less than a quarter of it apart meet head-on). `flat`: the +/// face is a plane, whose shading does not depend on its interior: a butterfly's square sides +/// take half the facing arc's count there (each wing steps 2:1), instead of the whole count. `None` when the face has no hole, or a hole has no room +/// for a collar (it touches another feature). +pub(crate) fn segments(uvs: &[[f64; 2]], loops: &[std::ops::Range], scale: [f64; 2], seed: f64, flat: bool) -> Option> { + let snap = 0.25 * seed; + let metric = |q: [f64; 2]| [q[0] * scale[0], q[1] * scale[1]]; + let unmetric = |q: [f64; 2]| [q[0] / scale[0], q[1] / scale[1]]; + let area = |r: &std::ops::Range| -> f64 { + let m = r.len(); + (0..m) + .filter_map(|k| { + let (a, b) = (metric(*uvs.get(r.start + k)?), metric(*uvs.get(r.start + (k + 1) % m)?)); + Some(0.5 * (a[0] * b[1] - b[0] * a[1])) + }) + .sum() + }; + // The outer loop encloses the most area; the others are holes. + let outer = loops.iter().enumerate().max_by(|a, b| area(a.1).abs().total_cmp(&area(b.1).abs()))?.0; + let holes: Vec<&std::ops::Range> = loops.iter().enumerate().filter(|(i, _)| *i != outer).map(|(_, r)| r).collect(); + let outer = loops.get(outer)?; + let bbox = |r: &std::ops::Range| -> Option<([f64; 2], [f64; 2])> { + let (mut lo, mut hi) = ([f64::MAX; 2], [f64::MIN; 2]); + for k in r.clone() { + let q = metric(*uvs.get(k)?); + for d in 0..2 { + lo[d] = lo[d].min(q[d]); + hi[d] = hi[d].max(q[d]); + } + } + Some((lo, hi)) + }; + // The node of loop `r` furthest out towards each corner of its box, in the box's own + // proportions (its 45° points on a round loop), corner by corner from (min, min). + let supports = |r: &std::ops::Range| -> Option<[usize; 4]> { + let (lo, hi) = bbox(r)?; + let (w, h) = (hi[0] - lo[0], hi[1] - lo[1]); + if !(w > 0.0 && h > 0.0) { + return None; + } + let mut out = [0usize; 4]; + for (slot, dir) in out.iter_mut().zip([[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]]) { + let reach = |i: usize| { + uvs.get(i).map_or(f64::MIN, |q| { + let q = metric(*q); + dir[0] * (q[0] - lo[0]) / w + dir[1] * (q[1] - lo[1]) / h + }) + }; + *slot = (r.start..r.end).max_by(|&x, &y| reach(x).total_cmp(&reach(y)))?; + } + let mut distinct = out.to_vec(); + distinct.sort_unstable(); + distinct.dedup(); + (distinct.len() == 4).then_some(out) + }; + // The arc of loop `r` (in nodes) between its support points `a` and `b`, the way that passes + // neither of the other two. + let side_count = |n: usize| if flat { n.div_ceil(2).max(1) } else { n.max(1) }; + let arc = |r: &std::ops::Range, p: &[usize; 4], k: usize| -> Option { + let m = r.len(); + let (a, b) = (*p.get(k)? - r.start, *p.get((k + 1) % 4)? - r.start); + let fwd = (b + m - a) % m; + let others: Vec = (0..4).filter(|&j| j != k && j != (k + 1) % 4).filter_map(|j| p.get(j).map(|x| (x - r.start + m - a) % m)).collect(); + Some(if others.iter().all(|&o| o > fwd) { fwd } else { m - fwd }) + }; + // A face of one round loop (no hole, no concave peak; a disc, an ellipse): the butterfly. + // A square inside it, its corners halfway from the loop's centre to the loop's four 45° + // points, each joined to its point; each square side in as many segments as the arc facing + // it. Five four-sided patches: the square and four wings. + if holes.is_empty() { + let alone: Vec<[f64; 2]> = outer.clone().filter_map(|k| uvs.get(k).copied()).collect(); + if !motorcycle::traces(&alone, &[std::ops::Range { start: 0, end: alone.len() }], scale, snap).is_empty() { + return None; + } + let p = supports(outer)?; + let (lo, hi) = bbox(outer)?; + let centre = [0.5 * (lo[0] + hi[0]), 0.5 * (lo[1] + hi[1])]; + let mut corners = [[0.0; 2]; 4]; + for (c, &i) in corners.iter_mut().zip(&p) { + let q = metric(*uvs.get(i)?); + *c = unmetric([centre[0] + 0.5 * (q[0] - centre[0]), centre[1] + 0.5 * (q[1] - centre[1])]); + } + let mut segs = Vec::with_capacity(8); + for k in 0..4 { + segs.push(Segment { + a: End::Point(*corners.get(k)?), + b: End::Point(*corners.get((k + 1) % 4)?), + count: Some(side_count(arc(outer, &p, k)?)), + }); + } + for (k, &i) in p.iter().enumerate() { + segs.push(Segment { a: End::Point(*corners.get(k)?), b: End::Node(u32::try_from(i).ok()?), count: None }); + } + return Some(segs); + } + // A square needs a strip at least a seed wide between it and the outer loop on every side + // (and a seed of collar inside it): a hole with less than two seeds to the outer loop on any + // side has no room for one, so the outer loop is its square, and the collar cuts join the + // hole's 45° points to the outer loop's. + let (olo, ohi) = bbox(outer)?; + for hole in &holes { + let (lo, hi) = bbox(hole)?; + let gaps = [lo[0] - olo[0], lo[1] - olo[1], ohi[0] - hi[0], ohi[1] - hi[1]]; + if holes.len() == 1 && gaps.iter().any(|g| *g < 2.0 * seed) { + let (inner, out) = (supports(hole)?, supports(outer)?); + return inner + .iter() + .zip(out) + .map(|(&a, b)| Some(Segment { a: End::Node(u32::try_from(a).ok()?), b: End::Node(u32::try_from(b).ok()?), count: None })) + .collect(); + } + } + let boxes: Vec<([f64; 2], [f64; 2])> = holes.iter().map(|r| bbox(r)).collect::>()?; + // Distances in the metric: point to box, point to segment, box to box. + let to_box = |p: [f64; 2], (lo, hi): ([f64; 2], [f64; 2])| { + let dx = (lo[0] - p[0]).max(0.0).max(p[0] - hi[0]); + let dy = (lo[1] - p[1]).max(0.0).max(p[1] - hi[1]); + dx.hypot(dy) + }; + let to_segment = |p: [f64; 2], a: [f64; 2], b: [f64; 2]| { + let (ab, ap) = ([b[0] - a[0], b[1] - a[1]], [p[0] - a[0], p[1] - a[1]]); + let len2 = ab[0] * ab[0] + ab[1] * ab[1]; + let t = if len2 > 0.0 { ((ap[0] * ab[0] + ap[1] * ab[1]) / len2).clamp(0.0, 1.0) } else { 0.0 }; + (p[0] - a[0] - t * ab[0]).hypot(p[1] - a[1] - t * ab[1]) + }; + let corners_of = |(lo, hi): ([f64; 2], [f64; 2])| [[lo[0], lo[1]], [hi[0], lo[1]], [hi[0], hi[1]], [lo[0], hi[1]]]; + // Each square: the hole's box grown by half its clearance to the outer loop and a quarter of + // its gap to another hole's box (both squares grow into that gap), so squares never touch. + let m = outer.len(); + let mut squares: Vec<([f64; 2], [f64; 2])> = Vec::with_capacity(boxes.len()); + for (i, &bx) in boxes.iter().enumerate() { + let mut clear = f64::INFINITY; + for k in 0..m { + let (a, b) = (metric(*uvs.get(outer.start + k)?), metric(*uvs.get(outer.start + (k + 1) % m)?)); + clear = clear.min(to_box(a, bx)); + for c in corners_of(bx) { + clear = clear.min(to_segment(c, a, b)); + } + } + for (j, &other) in boxes.iter().enumerate() { + if j != i { + // Box to box: the gap along each axis. + let gx = (other.0[0] - bx.1[0]).max(bx.0[0] - other.1[0]).max(0.0); + let gy = (other.0[1] - bx.1[1]).max(bx.0[1] - other.1[1]).max(0.0); + // Both holes' squares grow into this gap: each takes half of it. + clear = clear.min(0.5 * gx.hypot(gy)); + } + } + let margin = 0.5 * clear; + if !(margin.is_finite() && margin > 0.0) { + return None; + } + squares.push(([bx.0[0] - margin, bx.0[1] - margin], [bx.1[0] + margin, bx.1[1] + margin])); + } + let mut segs: Vec = Vec::new(); + // Each hole's four 45° points (towards its square's corners, in order). + let picks: Vec<[usize; 4]> = holes.iter().map(|r| supports(r)).collect::>()?; + // Square sides (four per hole, corner k to corner k + 1), each in as many segments as the + // hole's arc facing it (between its 45° points k and k + 1, the way that passes neither + // other point): the wings of the butterfly then have matching opposite sides. + for ((&sq, hole), p) in squares.iter().zip(&holes).zip(&picks) { + let c = corners_of(sq).map(unmetric); + for k in 0..4 { + segs.push(Segment { a: End::Point(*c.get(k)?), b: End::Point(*c.get((k + 1) % 4)?), count: Some(side_count(arc(hole, p, k)?)) }); + } + } + for (p, &sq) in picks.iter().zip(&squares) { + let c = corners_of(sq).map(unmetric); + for (k, &i) in p.iter().enumerate() { + segs.push(Segment { a: End::Point(*c.get(k)?), b: End::Node(u32::try_from(i).ok()?), count: None }); + } + } + // The motorcycles, with the squares as the obstacles: the outer loop, then each square as a + // loop of its four corners. + let mut domain: Vec<[f64; 2]> = outer.clone().filter_map(|k| uvs.get(k).copied()).collect(); + let mut dloops = vec![std::ops::Range { start: 0, end: domain.len() }]; + for &sq in &squares { + let start = domain.len(); + domain.extend(corners_of(sq).map(unmetric)); + dloops.push(start..start + 4); + } + let traces = motorcycle::traces(&domain, &dloops, scale, snap); + let base = segs.len(); + let on_outer = |i: usize| i < outer.len(); + for (ti, t) in traces.iter().enumerate() { + let start = if on_outer(t.start) { End::Node(u32::try_from(outer.start + t.start).ok()?) } else { End::Point(*domain.get(t.start)?) }; + let end = match t.end { + Stop::Boundary(a, b) if on_outer(a) && on_outer(b) => { + // The nearer end of the boundary segment it reached (a cut never adds a node to + // the boundary). + let (pa, pb) = (metric(*domain.get(a)?), metric(*domain.get(b)?)); + let q = metric(t.b); + let k = if (pa[0] - q[0]).hypot(pa[1] - q[1]) <= (pb[0] - q[0]).hypot(pb[1] - q[1]) { a } else { b }; + End::Node(u32::try_from(outer.start + k).ok()?) + } + Stop::Boundary(a, b) => { + // A side of a square: corners `a` and `b` of square `s`. + let (s, la, lb) = ((a - outer.len()) / 4, (a - outer.len()) % 4, (b - outer.len()) % 4); + let k = if (la + 1) % 4 == lb { la } else { lb }; + End::On(4 * s + k, t.b) + } + Stop::Trace(c) => End::On(base + c, t.b), + Stop::HeadOn(j) => { + let o = traces.get(j)?; + End::Point([0.5 * (t.b[0] + o.b[0]), 0.5 * (t.b[1] + o.b[1])]) + } + }; + let _ = ti; + segs.push(Segment { a: start, b: end, count: None }); + } + Some(segs) +} diff --git a/crates/kernel/src/tess/coons.rs b/crates/kernel/src/tess/coons.rs index 6f79d518..e803f7dc 100644 --- a/crates/kernel/src/tess/coons.rs +++ b/crates/kernel/src/tess/coons.rs @@ -14,11 +14,20 @@ const CORNER_TURN: f64 = 0.6; struct Side { nodes: Vec, at: Vec, + /// Its length in 3D (mm). + len: f64, } /// The four corners of a single boundary loop, by how sharply the loop turns there in (u, v) /// scaled by `scale` (one unit per tolerance cell); `None` unless there are exactly four. fn corners(uv: &[[f64; 2]], scale: [f64; 2]) -> Option<[usize; 4]> { + let [a, b, c, d] = candidates(uv, scale)?[..] else { return None }; + Some([a, b, c, d]) +} + +/// Every node where the loop turns more than [`CORNER_TURN`] (or meets a pole copy), in loop +/// order. +fn candidates(uv: &[[f64; 2]], scale: [f64; 2]) -> Option> { let n = uv.len(); if n < 4 { return None; @@ -39,8 +48,60 @@ fn corners(uv: &[[f64; 2]], scale: [f64; 2]) -> Option<[usize; 4]> { found.push(i); } } - let [a, b, c, d] = found[..] else { return None }; - Some([a, b, c, d]) + Some(found) +} + +/// Most segments a patch's own level gives one direction. +const MAX_LEVEL: usize = 4096; + +/// Most corner candidates [`coons_best`] chooses four from (8 choose 4 = 70 blends at most). +const MAX_CANDIDATES: usize = 8; + +/// Like [`coons`] for a loop with more than four corners (a patch of a trimmed face with a +/// corner cut off, the wing of a butterfly around a hole): every choice of four of its corners +/// is blended, and the one whose smallest triangle angle (in the metric `scale`) is largest +/// wins. Its other corners lie on sides. `None` if no choice gives a blend that does not fold. +pub(crate) fn coons_best(uv: &[[f64; 2]], fixed: &[Vec3], scale: [f64; 2], max_points: usize, h: f64) -> Option<(Vec<[f64; 2]>, Vec<[u32; 3]>)> { + let c = candidates(uv, scale)?; + if c.len() < 4 || c.len() > MAX_CANDIDATES { + return None; + } + let mut best: Option<(f64, (Vec<[f64; 2]>, Vec<[u32; 3]>))> = None; + let k = c.len(); + for a in 0..k { + for b in a + 1..k { + for d in b + 1..k { + for e in d + 1..k { + let pick = [*c.get(a)?, *c.get(b)?, *c.get(d)?, *c.get(e)?]; + let Some(r) = coons_at(uv, fixed, pick, max_points, h, None) else { continue }; + let q = min_angle(uv, &r, scale); + if best.as_ref().is_none_or(|(bq, _)| q > *bq) { + best = Some((q, r)); + } + } + } + } + } + best.map(|(_, r)| r) +} + +/// Smallest corner angle (radians) of the triangles, in the metric. +fn min_angle(uv: &[[f64; 2]], r: &(Vec<[f64; 2]>, Vec<[u32; 3]>), scale: [f64; 2]) -> f64 { + let at = |i: u32| { + let i = i as usize; + let q = if i < uv.len() { uv.get(i) } else { r.0.get(i - uv.len()) }; + q.map(|q| [q[0] * scale[0], q[1] * scale[1]]) + }; + let mut worst = std::f64::consts::PI; + for t in &r.1 { + let (Some(a), Some(b), Some(c)) = (at(t[0]), at(t[1]), at(t[2])) else { continue }; + for (p, q, s) in [(a, b, c), (b, c, a), (c, a, b)] { + let (u, v) = ([q[0] - p[0], q[1] - p[1]], [s[0] - p[0], s[1] - p[1]]); + let ang = (u[0] * v[1] - u[1] * v[0]).abs().atan2(u[0] * v[0] + u[1] * v[1]); + worst = worst.min(ang); + } + } + worst } fn side(fixed: &[Vec3], n: usize, from: usize, to: usize) -> Option { @@ -58,9 +119,9 @@ fn side(fixed: &[Vec3], n: usize, from: usize, to: usize) -> Option { if !(total > 0.0) { // A pole side: evenly along it. let m = (nodes.len() - 1).max(1) as f64; - return Some(Side { at: (0..nodes.len()).map(|k| k as f64 / m).collect(), nodes }); + return Some(Side { at: (0..nodes.len()).map(|k| k as f64 / m).collect(), nodes, len: 0.0 }); } - Some(Side { nodes, at: at.into_iter().map(|x| x / total).collect() }) + Some(Side { nodes, at: at.into_iter().map(|x| x / total).collect(), len: total }) } /// The point of a side at position `s` (0 to 1), linear between its nodes. @@ -72,9 +133,12 @@ fn along(side: &Side, uv: &[[f64; 2]], s: f64) -> Option<[f64; 2]> { Some([p[0] + (q[0] - p[0]) * f, p[1] + (q[1] - p[1]) * f]) } -/// Interior line positions for a pair of opposite sides: the sides' own when they have the -/// same nodes (so the grid meets both node for node), else evenly as many as the denser side. -fn lines(a: &Side, b: &Side) -> Vec { +/// Interior line positions for a pair of opposite sides: the sides' own when they have as many +/// nodes as each other (so the grid meets both node for node); else evenly at the inner level a +/// GPU tessellator takes, the larger of the two sides' (outer) levels, and never coarser than the +/// size rules ask (the sides' mean length over the target side `h`). The denser side then meets +/// the grid node for node and the sparser one steps down to it through its ring. +fn lines(a: &Side, b: &Side, h: f64) -> Vec { let inner = |s: &Side| s.at.get(1..s.at.len().saturating_sub(1)).map(<[f64]>::to_vec).unwrap_or_default(); if a.at.len() == b.at.len() { let l: Vec = inner(a).iter().zip(inner(b)).map(|(x, y)| 0.5 * (x + y)).collect(); @@ -82,23 +146,43 @@ fn lines(a: &Side, b: &Side) -> Vec { return l; } } - let n = a.at.len().max(b.at.len()).saturating_sub(1).max(2); + let size = if h > 0.0 { (0.5 * (a.len + b.len) / h).ceil() } else { 1.0 }; + let size = if size.is_finite() { size as usize } else { 1 }; + let outer = a.at.len().max(b.at.len()).saturating_sub(1); + let n = outer.max(size).clamp(2, MAX_LEVEL); (1..n).map(|i| i as f64 / n as f64).collect() } /// Triangulate the loop `uv` (with its fixed 3D nodes) as a Coons patch: boundary vertices /// keep their indices, the interior grid follows. Triangles counter-clockwise in (u, v). -/// `fine(a, b)`: whether the surface between two parameter points needs no point between them -/// (the grid lines, taken from the sides, are split inside the face where it does not). +/// +/// `fine(a, b)`: whether the surface between two parameter points needs no point between them. +/// Given for a whole face, whose sides are its edges' nodes: where its grid converges (towards a +/// pole) the interior can need more lines than the sides carry, and the lines are split there. +/// A motorcycle patch passes `None`: its sides are cuts seeded at the face's global seed. pub(crate) fn coons( uv: &[[f64; 2]], fixed: &[Vec3], scale: [f64; 2], max_points: usize, - fine: &dyn Fn([f64; 2], [f64; 2]) -> bool, + h: f64, + fine: Option<&dyn Fn([f64; 2], [f64; 2]) -> bool>, +) -> Option<(Vec<[f64; 2]>, Vec<[u32; 3]>)> { + let pick = corners(uv, scale)?; + coons_at(uv, fixed, pick, max_points, h, fine) +} + +/// [`coons`] with the four corners given (loop indices, in loop order). +pub(crate) fn coons_at( + uv: &[[f64; 2]], + fixed: &[Vec3], + pick: [usize; 4], + max_points: usize, + h: f64, + fine: Option<&dyn Fn([f64; 2], [f64; 2]) -> bool>, ) -> Option<(Vec<[f64; 2]>, Vec<[u32; 3]>)> { let n = uv.len(); - let [c0, c1, c2, c3] = corners(uv, scale)?; + let [c0, c1, c2, c3] = pick; // Sides in loop order; the bottom (0) and right (1) run with the loop, the top (2) and // left (3) against it, so s and t both run from corner c0. let bottom = side(fixed, n, c0, c1)?; @@ -116,8 +200,16 @@ pub(crate) fn coons( }; Some([blend(0), blend(1)]) }; - let ss = refine(lines(&bottom, &top), max_points, &|a, b, x| Some(fine(at(a, x)?, at(b, x)?))); - let ts = refine(lines(&left, &right), max_points, &|a, b, x| Some(fine(at(x, a)?, at(x, b)?))); + // The lines come from the sides' seeds (see [`lines`]): every side is seeded at the face's + // global seed, which the size rules set for the whole surface, so the interior needs no finer + // spacing than its sides. + let (ss, ts) = (lines(&bottom, &top, h), lines(&left, &right, h)); + let (ss, ts) = match fine { + Some(fine) => { + (refine(ss, max_points, &|a, b, x| Some(fine(at(a, x)?, at(b, x)?))), refine(ts, max_points, &|a, b, x| Some(fine(at(x, a)?, at(x, b)?)))) + } + None => (ss, ts), + }; if ss.len().checked_mul(ts.len())? > max_points { return None; } @@ -141,10 +233,46 @@ pub(crate) fn coons( // the bottom and right strips come out counter-clockwise, the top and left clockwise. let row = |j: usize| (0..ns).map(|i| g(i, j)).collect::>(); let col = |i: usize| (0..nt).map(|j| g(i, j)).collect::>(); - zip(&bottom, &row(0), &ss, false, &mut tris)?; - zip(&right, &col(ns - 1), &ts, false, &mut tris)?; - zip(&top, &row(nt - 1), &ss, true, &mut tris)?; - zip(&left, &col(0), &ts, true, &mut tris)?; + // Where a side's segments and its grid row's differ, a transition row steps between them + // halfway to the first grid line (its segments the two counts' geometric mean), each zipped + // to it, with one triangle at each end closing it against the corner. + let (s_first, s_last) = (*ss.first()?, *ss.last()?); + let (t_first, t_last) = (*ts.first()?, *ts.last()?); + // Each side: its grid row, the row's positions, the orientation, where the transition row + // runs, and the point at a position along it. The row lies halfway between the side and the + // grid row; its ends lie halfway between the diagonal from the corner to the grid row's end + // and that end, so they are strictly inside the side's own ring. + #[allow(clippy::type_complexity)] + let sides: [(&Side, Vec, &[f64], bool, (f64, f64), &dyn Fn(f64) -> [f64; 2]); 4] = [ + (&bottom, row(0), &ss, false, (0.75 * s_first, 1.0 - 0.75 * (1.0 - s_last)), &|x| [x, 0.5 * t_first]), + (&right, col(ns - 1), &ts, false, (0.75 * t_first, 1.0 - 0.75 * (1.0 - t_last)), &|x| [0.5 * (1.0 + s_last), x]), + (&top, row(nt - 1), &ss, true, (0.75 * s_first, 1.0 - 0.75 * (1.0 - s_last)), &|x| [x, 0.5 * (1.0 + t_last)]), + (&left, col(0), &ts, true, (0.75 * t_first, 1.0 - 0.75 * (1.0 - t_last)), &|x| [0.5 * s_first, x]), + ]; + for (side, inner, inner_at, flip, (x0, x1), place) in sides { + let (outer_segs, inner_segs) = (side.nodes.len().saturating_sub(1), inner.len() + 1); + let mean = ((outer_segs * inner_segs) as f64).sqrt().round() as usize; + if mean == outer_segs || mean == inner_segs || mean < 2 { + zip(side, &inner, inner_at, flip, &mut tris)?; + continue; + } + let k = mean - 1; + let at_t: Vec = (0..k).map(|j| if k > 1 { x0 + (x1 - x0) * j as f64 / (k - 1) as f64 } else { 0.5 * (x0 + x1) }).collect(); + let mut mid: Vec = Vec::with_capacity(k); + for &x in &at_t { + let [sx, tx] = place(x); + mid.push(u32::try_from(out.len()).ok()?); + out.push(at(sx, tx)?); + } + zip(side, &mid, &at_t, flip, &mut tris)?; + zip(&Side { nodes: mid.clone(), at: at_t.clone(), len: 0.0 }, &inner, inner_at, flip, &mut tris)?; + // The two end triangles: corner, transition end, grid end, turned like the zips' (in + // (s, t) the first runs with them and the last against). Never turned to fit the blend: + // a blend that folds there leaves one inverted, and the patch is refused below. + for t in [[*side.nodes.first()?, *mid.first()?, *inner.first()?], [*side.nodes.last()?, *inner.last()?, *mid.last()?]] { + tris.push(if flip { [t[0], t[2], t[1]] } else { t }); + } + } // The (s, t) square maps to (u, v) keeping or reversing orientation, with the loop: all // triangles must share one sign, which is then made counter-clockwise. let area = |t: &[u32; 3]| -> Option { @@ -188,7 +316,7 @@ fn refine(lines: Vec, max: usize, fine: &dyn Fn(f64, f64, f64) -> Option Side { - Side { nodes: s.nodes.into_iter().rev().collect(), at: s.at.into_iter().rev().map(|x| 1.0 - x).collect() } + Side { nodes: s.nodes.into_iter().rev().collect(), at: s.at.into_iter().rev().map(|x| 1.0 - x).collect(), len: s.len } } /// Join a side (its nodes at positions `side.at`) to a grid row at positions `inner_at`. @@ -251,7 +379,7 @@ mod tests { b.reverse(); } let fixed: Vec = b.iter().map(|p| Vec3::new(p[0], p[1], 0.0)).collect(); - let Some((extra, tris)) = coons(&b, &fixed, [1.0, 1.0], 10_000, &|_, _| true) else { panic!("meshes (k {k})") }; + let Some((extra, tris)) = coons(&b, &fixed, [1.0, 1.0], 10_000, 1.0, None) else { panic!("meshes (k {k})") }; let mut uv = b.clone(); uv.extend(extra); assert!((area(&uv, &tris) - 16.0).abs() < 1e-9, "k {k} reverse {reverse}: {}", area(&uv, &tris)); @@ -259,10 +387,146 @@ mod tests { } } + /// Opposite sides with different node counts (20 against 4): the interior grid takes the + /// larger side's level (a tessellator's inner level), so the dense side meets it node for node + /// and the sparse side steps down to it through a transition row. Covered exactly, nothing + /// missing or doubled. + #[test] + fn the_interior_level_is_the_larger_sides_and_the_sparser_side_steps_down() { + // A 4 × 4 square: bottom in 20 segments, the other sides in 4. + let mut b: Vec<[f64; 2]> = (0..20).map(|k| [k as f64 * 0.2, 0.0]).collect(); + b.extend((0..4).map(|k| [4.0, k as f64])); + b.extend((0..4).map(|k| [4.0 - k as f64, 4.0])); + b.extend((0..4).map(|k| [0.0, 4.0 - k as f64])); + let fixed: Vec = b.iter().map(|p| Vec3::new(p[0], p[1], 0.0)).collect(); + let Some((extra, tris)) = coons(&b, &fixed, [1.0, 1.0], 10_000, 1.0, None) else { panic!("meshes") }; + let mut uv = b.clone(); + uv.extend(extra.iter().copied()); + assert!((area(&uv, &tris) - 16.0).abs() < 1e-9); + // Level 20 along s (the bottom's), 4 along t: 19 × 3 interior points; plus the transition + // row by the 4-segment top (√(4 · 20) ≈ 9 segments: 8 points). + assert_eq!(extra.len(), 19 * 3 + 8, "{}", extra.len()); + // Every inner side shared by two triangles, every boundary side used once: nothing + // missing or doubled at the transition's ends. + let mut sides: std::collections::HashMap<(u32, u32), i32> = std::collections::HashMap::new(); + for t in &tris { + for e in 0..3 { + *sides.entry((t[e], t[(e + 1) % 3])).or_default() += 1; + } + } + let n = b.len() as u32; + for (&(a, c), &k) in &sides { + assert_eq!(k, 1, "side {a}-{c} once each way"); + let boundary = a < n && c < n && (c == (a + 1) % n); + assert!(boundary || sides.contains_key(&(c, a)), "inner side {a}-{c} shared"); + } + } + #[test] fn five_corners_is_not_a_coons_patch() { let b = [[0.0, 0.0], [1.0, 0.0], [1.0, 0.5], [0.5, 1.0], [0.0, 1.0]]; let fixed: Vec = b.iter().map(|p| Vec3::new(p[0], p[1], 0.0)).collect(); - assert!(coons(&b, &fixed, [1.0, 1.0], 100, &|_, _| true).is_none()); + assert!(coons(&b, &fixed, [1.0, 1.0], 100, 1.0, None).is_none()); + } + + /// A planar strip with a fillet at one end and a chamfer at the other, in mm (a face of + /// robustness seed 33): its corners make the one-segment chamfer a side facing the long one, + /// and at some target sides the blend folds beside the chamfer. The transition row's end + /// triangles are turned like its zips, never to fit the blend, so a fold leaves a triangle + /// inverted and the blend is refused; any blend given covers the loop exactly once. + #[test] + #[allow(clippy::approx_constant)] // measured coordinates, one of them near ln 10 + fn a_blend_folding_beside_a_short_side_is_refused_not_overlapped() { + let b: Vec<[f64; 2]> = vec![ + [24.0000, 3.0150], + [23.0000, 4.4292], + [16.6667, 4.4292], + [10.3333, 4.4292], + [4.0000, 4.4292], + [4.0000, 3.0150], + [4.0002, 2.9836], + [4.0010, 2.9520], + [4.0023, 2.9201], + [4.0040, 2.8880], + [4.0064, 2.8556], + [4.0092, 2.8231], + [4.0127, 2.7905], + [4.0167, 2.7577], + [4.0213, 2.7248], + [4.0265, 2.6919], + [4.0322, 2.6589], + [4.0386, 2.6259], + [4.0456, 2.5929], + [4.0532, 2.5600], + [4.0614, 2.5272], + [4.0702, 2.4945], + [4.0797, 2.4619], + [4.0897, 2.4295], + [4.1004, 2.3974], + [4.1117, 2.3655], + [4.1236, 2.3339], + [4.1362, 2.3026], + [4.1493, 2.2717], + [4.1630, 2.2412], + [4.1773, 2.2111], + [4.1922, 2.1814], + [4.2076, 2.1523], + [4.2236, 2.1237], + [4.2402, 2.0956], + [4.2573, 2.0681], + [4.2748, 2.0412], + [4.2929, 2.0150], + [4.3114, 1.9895], + [4.3304, 1.9646], + [4.3499, 1.9405], + [4.3697, 1.9171], + [4.3899, 1.8945], + [4.4105, 1.8726], + [4.4315, 1.8516], + [4.4528, 1.8313], + [4.4744, 1.8119], + [4.4962, 1.7934], + [4.5184, 1.7757], + [4.5407, 1.7588], + [4.5633, 1.7428], + [4.5860, 1.7277], + [4.6089, 1.7135], + [4.6319, 1.7001], + [4.6550, 1.6876], + [4.6783, 1.6760], + [4.7015, 1.6653], + [4.7248, 1.6554], + [4.7482, 1.6464], + [4.7715, 1.6382], + [4.7948, 1.6309], + [4.8180, 1.6244], + [4.8412, 1.6188], + [4.8643, 1.6139], + [4.8873, 1.6098], + [4.9101, 1.6065], + [4.9329, 1.6040], + [4.9554, 1.6022], + [4.9778, 1.6012], + [5.0000, 1.6008], + [11.3333, 1.6008], + [17.6667, 1.6008], + [24.0000, 1.6008], + ]; + let fixed: Vec = b.iter().map(|p| Vec3::new(p[0], p[1], 0.0)).collect(); + let n = b.len(); + let loop_area: f64 = (0..n).map(|k| 0.5 * (b[k][0] * b[(k + 1) % n][1] - b[(k + 1) % n][0] * b[k][1])).sum(); + for h in [0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 8.0] { + let Some((extra, tris)) = coons(&b, &fixed, [1.0, 1.0], 100_000, h, None) else { continue }; + let mut uv = b.clone(); + uv.extend(extra); + assert!((area(&uv, &tris) - loop_area).abs() < 1e-9 * loop_area, "h {h}: {} for {loop_area}", area(&uv, &tris)); + let mut sides: std::collections::HashMap<(u32, u32), usize> = std::collections::HashMap::new(); + for t in &tris { + for k in 0..3 { + *sides.entry((t[k], t[(k + 1) % 3])).or_default() += 1; + } + } + assert!(sides.values().all(|&k| k == 1), "h {h}: a side used twice the same way"); + } } } diff --git a/crates/kernel/src/tess/cutgrid.rs b/crates/kernel/src/tess/cutgrid.rs new file mode 100644 index 00000000..2c53dc8a --- /dev/null +++ b/crates/kernel/src/tess/cutgrid.rs @@ -0,0 +1,518 @@ +//! A trimmed face as an iso-line T-mesh: full-length grid lines from the size rules plus the +//! face's motorcycle traces (which run only as far as they reach), whose rectangles that lie +//! well inside the face are meshed directly. The band between those rectangles (a staircase) +//! and the face's boundary nodes is triangulated with both kept as constraints and no point +//! added. The boundary stays exactly the shared edge nodes, so the face meets its neighbours +//! exactly, and only the band next to the boundary is irregular (cut-cell meshing). + +use super::cdt; +use super::motorcycle::Trace; + +/// A cell is kept when the boundary stays at least this fraction of its smaller side away from +/// it (in the metric): closer, the band triangle between them would be a sliver. +const CLEARANCE: f64 = 0.35; +/// Two lines closer than this fraction of a grid cell would leave a strip of slivers between +/// them: a grid line that close to a trace moves onto the trace (the trace runs through its +/// peak, which is the point of it), and a trace that close to an earlier one joins it. +const SNAP: f64 = 0.25; + +/// One direction's lines: positions (increasing) and whether each runs the whole face. +struct Lines { + at: Vec, + full: Vec, +} + +/// The full lines `grid` plus the traces' lines in this direction (see [`SNAP`]; the first and +/// last grid lines, the face's parameter box, never move). Returns the lines and, per trace +/// (in order), its line's index. +fn lines(grid: &[f64], traces: &[f64]) -> Option<(Lines, Vec)> { + let mut at: Vec = grid.to_vec(); + let mut full = vec![true; at.len()]; + let mut place = Vec::with_capacity(traces.len()); + for &c in traces { + // The full-grid cell around c, and every line inside it. + let k = grid.partition_point(|g| *g < c); + let (lo, hi) = (*grid.get(k.saturating_sub(1))?, *grid.get(k.min(grid.len().saturating_sub(1)))?); + let w = (hi - lo).max(f64::MIN_POSITIVE); + let last = grid.len().saturating_sub(1); + let near = at.iter().enumerate().filter(|(_, x)| (**x - c).abs() < SNAP * w).min_by(|a, b| (a.1 - c).abs().total_cmp(&(b.1 - c).abs())); + match near { + // An earlier trace, or a box edge: join it. + Some((i, _)) if !full.get(i).copied().unwrap_or(false) || i == 0 || i == last => place.push(i), + // An inner grid line: it moves onto the trace. + Some((i, _)) => { + *at.get_mut(i)? = c; + place.push(i); + } + None => { + at.push(c); + full.push(false); + place.push(at.len() - 1); + } + } + } + // Sort, carrying the trace indices along. + let mut order: Vec = (0..at.len()).collect(); + order.sort_by(|a, b| at.get(*a).copied().unwrap_or(0.0).total_cmp(&at.get(*b).copied().unwrap_or(0.0))); + let mut rank = vec![0; at.len()]; + for (r, &i) in order.iter().enumerate() { + *rank.get_mut(i)? = r; + } + let lines = + Lines { at: order.iter().filter_map(|&i| at.get(i).copied()).collect(), full: order.iter().filter_map(|&i| full.get(i).copied()).collect() }; + let place = place.iter().filter_map(|&i| rank.get(i).copied()).collect::>(); + (place.len() == traces.len()).then_some((lines, place)) +} + +/// Mesh the region bounded by `loops` (ranges of `boundary`, each closed; even–odd, so inner +/// loops are holes) with the full grid lines `us` × `vs` (increasing, spanning the region's +/// parameter box) and the motorcycle `traces`, measured in the metric `scale`. Boundary +/// vertices keep their indices; grid points follow. Triangles counter-clockwise in (u, v). +/// `None` when no cell is kept or the band cannot be triangulated (the face then goes to +/// Delaunay). +pub(crate) fn cut_grid( + boundary: &[[f64; 2]], + loops: &[std::ops::Range], + us: &[f64], + vs: &[f64], + traces: &[Trace], + scale: [f64; 2], +) -> Option<(Vec<[f64; 2]>, Vec<[u32; 3]>)> { + if us.len() < 2 || vs.len() < 2 || !scale.iter().all(|s| s.is_finite() && *s > 0.0) { + return None; + } + // Traces on u = c (axis 0) are lines in u; on v = c, lines in v. + let tu: Vec<&Trace> = traces.iter().filter(|t| t.axis == 0).collect(); + let tv: Vec<&Trace> = traces.iter().filter(|t| t.axis == 1).collect(); + let (lu, pu) = lines(us, &tu.iter().map(|t| t.at).collect::>())?; + let (lv, pv) = lines(vs, &tv.iter().map(|t| t.at).collect::>())?; + let (nu, nv) = (lu.at.len(), lv.at.len()); + let (cu, cv) = (nu - 1, nv - 1); + let cells = cu.checked_mul(cv)?; + if cells > super::grid::MAX_FACE_GRID { + return None; + } + // Which sides of the fine cells are lines: the vertical side at column i in row j, the + // horizontal side at row j in column i. A trace covers a side it runs at least half along. + let covers = |t: &Trace, a: f64, b: f64| (t.to.min(b) - t.from.max(a)) >= 0.5 * (b - a); + let mut vert = vec![false; nu * cv]; + for i in 0..nu { + let full = lu.full.get(i).copied().unwrap_or(false); + for j in 0..cv { + let (a, b) = (*lv.at.get(j)?, *lv.at.get(j + 1)?); + *vert.get_mut(j * nu + i)? = full || tu.iter().zip(&pu).any(|(t, &k)| k == i && covers(t, a, b)); + } + } + let mut horiz = vec![false; nv * cu]; + for j in 0..nv { + let full = lv.full.get(j).copied().unwrap_or(false); + for i in 0..cu { + let (a, b) = (*lu.at.get(i)?, *lu.at.get(i + 1)?); + *horiz.get_mut(j * cu + i)? = full || tv.iter().zip(&pv).any(|(t, &k)| k == j && covers(t, a, b)); + } + } + // Rectangles: fine cells joined across sides that are no line. A T-mesh's faces are + // rectangles; anything else (a trace ending in the open after snapping) falls back to + // every line running the whole face. + let rects = rectangles(cu, cv, &vert, &horiz).or_else(|| { + vert.iter_mut().for_each(|s| *s = true); + horiz.iter_mut().for_each(|s| *s = true); + rectangles(cu, cv, &vert, &horiz) + })?; + // In the metric, so clearances are lengths on the surface. + let gx: Vec = lu.at.iter().map(|u| u * scale[0]).collect(); + let gy: Vec = lv.at.iter().map(|v| v * scale[1]).collect(); + let segs: Vec<([f64; 2], [f64; 2])> = loops + .iter() + .flat_map(|r| { + let n = r.len(); + (0..n).filter_map(move |k| Some((*boundary.get(r.start + k)?, *boundary.get(r.start + (k + 1) % n)?))) + }) + .map(|(a, b)| ([a[0] * scale[0], a[1] * scale[1]], [b[0] * scale[0], b[1] * scale[1]])) + .collect(); + let keep = clear_cells(&gx, &gy, &segs, &rects)?; + if !keep.iter().any(|k| *k) { + return None; + } + let kept = |i: isize, j: isize| -> bool { + i >= 0 + && j >= 0 + && (i as usize) < cu + && (j as usize) < cv + && rects.of.get(j as usize * cu + i as usize).and_then(|r| keep.get(*r)).copied().unwrap_or(false) + }; + // Grid points on the kept rectangles' corners, numbered after the boundary. + let mut used = vec![false; nu * nv]; + for (r, b) in rects.boxes.iter().enumerate() { + if keep.get(r).copied().unwrap_or(false) { + for (i, j) in [(b[0], b[1]), (b[2], b[1]), (b[2], b[3]), (b[0], b[3])] { + *used.get_mut(j * nu + i)? = true; + } + } + } + let n = boundary.len(); + let mut id: Vec> = vec![None; nu * nv]; + let mut extra: Vec<[f64; 2]> = Vec::new(); + let mut node = |i: usize, j: usize, extra: &mut Vec<[f64; 2]>| -> Option { + let slot = id.get_mut(j * nu + i)?; + if let Some(k) = *slot { + return Some(k); + } + let k = u32::try_from(n + extra.len()).ok()?; + extra.push([*lu.at.get(i)?, *lv.at.get(j)?]); + *slot = Some(k); + Some(k) + }; + let mut tris: Vec<[u32; 3]> = Vec::new(); + // The staircase: pieces of kept rectangles' sides with no kept rectangle across, + // counter-clockwise around the kept region (so the band lies to their right). + let mut stairs: Vec<(u32, u32)> = Vec::new(); + for (r, b) in rects.boxes.iter().enumerate() { + if !keep.get(r).copied().unwrap_or(false) { + continue; + } + let [i0, j0, i1, j1] = *b; + // The perimeter counter-clockwise as grid points, each with the fine cell across the + // side that leaves it. + let mut ring: Vec<((usize, usize), (isize, isize))> = Vec::new(); + for i in i0..i1 { + ring.push(((i, j0), (i as isize, j0 as isize - 1))); + } + for j in j0..j1 { + ring.push(((i1, j), (i1 as isize, j as isize))); + } + for i in (i0 + 1..=i1).rev() { + ring.push(((i, j1), (i as isize - 1, j1 as isize))); + } + for j in (j0 + 1..=j1).rev() { + ring.push(((i0, j), (i0 as isize - 1, j as isize - 1))); + } + let corners = ring.iter().filter(|((i, j), _)| used.get(j * nu + i).copied().unwrap_or(false)).count(); + let mut poly: Vec = Vec::with_capacity(corners); + let mut across: Vec<(isize, isize)> = Vec::with_capacity(corners); + for &((i, j), c) in &ring { + if used.get(j * nu + i).copied().unwrap_or(false) { + poly.push(node(i, j, &mut extra)?); + across.push(c); + } + } + if poly.len() == 4 { + let [a, b, c, d] = [*poly.first()?, *poly.get(1)?, *poly.get(2)?, *poly.get(3)?]; + tris.push([a, b, c]); + tris.push([a, c, d]); + } else { + // Hanging nodes from neighbours: a fan from the rectangle's centre. + let centre = [0.5 * (lu.at.get(i0)? + lu.at.get(i1)?), 0.5 * (lv.at.get(j0)? + lv.at.get(j1)?)]; + let m = u32::try_from(n + extra.len()).ok()?; + extra.push(centre); + for k in 0..poly.len() { + tris.push([m, *poly.get(k)?, *poly.get((k + 1) % poly.len())?]); + } + } + for k in 0..poly.len() { + let (ci, cj) = *across.get(k)?; + if !kept(ci, cj) { + stairs.push((*poly.get(k)?, *poly.get((k + 1) % poly.len())?)); + } + } + } + let mut all = boundary.to_vec(); + all.extend_from_slice(&extra); + let band = cdt::band(&all, n, loops, &stairs, scale)?; + tris.extend(band); + // Every triangle counter-clockwise with some area: one that is not has folded. + for t in &tris { + let (a, b, c) = (all.get(t[0] as usize)?, all.get(t[1] as usize)?, all.get(t[2] as usize)?); + let area = (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]); + if !(area > 0.0) { + return None; + } + } + Some((extra, tris)) +} + +/// The rectangles of a T-mesh over `cu` × `cv` fine cells: per cell its rectangle, and per +/// rectangle its corners as grid indices `[i0, j0, i1, j1]`. `None` if a face is not a +/// rectangle. +struct Rects { + of: Vec, + boxes: Vec<[usize; 4]>, +} + +fn rectangles(cu: usize, cv: usize, vert: &[bool], horiz: &[bool]) -> Option { + let nu = cu + 1; + let mut of = vec![usize::MAX; cu * cv]; + let mut boxes = Vec::new(); + for j0 in 0..cv { + for i0 in 0..cu { + if *of.get(j0 * cu + i0)? != usize::MAX { + continue; + } + // Grow right while no line separates, then up while the whole row is open. + let mut i1 = i0 + 1; + while i1 < cu && !*vert.get(j0 * nu + i1)? { + i1 += 1; + } + let mut j1 = j0 + 1; + while j1 < cv && (i0..i1).all(|i| horiz.get(j1 * cu + i).is_some_and(|h| !*h)) { + j1 += 1; + } + let r = boxes.len(); + for j in j0..j1 { + for i in i0..i1 { + let slot = of.get_mut(j * cu + i)?; + if *slot != usize::MAX { + return None; + } + *slot = r; + } + } + boxes.push([i0, j0, i1, j1]); + } + } + // A rectangle: no line inside it, and lines all round it. + for b in &boxes { + let [i0, j0, i1, j1] = *b; + for j in j0..j1 { + for i in i0 + 1..i1 { + if *vert.get(j * nu + i)? { + return None; + } + } + } + for i in i0..i1 { + for j in j0 + 1..j1 { + if *horiz.get(j * cu + i)? { + return None; + } + } + } + for j in j0..j1 { + if (i0 > 0 && !*vert.get(j * nu + i0)?) || (i1 < cu && !*vert.get(j * nu + i1)?) { + return None; + } + } + for i in i0..i1 { + if (j0 > 0 && !*horiz.get(j0 * cu + i)?) || (j1 < cv && !*horiz.get(j1 * cu + i)?) { + return None; + } + } + } + Some(Rects { of, boxes }) +} + +/// Per rectangle: kept when every fine cell in it is inside the region and clear of the +/// boundary (see [`CLEARANCE`]). `gx`, `gy`: the lines in the metric; `segs`: the boundary. +fn clear_cells(gx: &[f64], gy: &[f64], segs: &[([f64; 2], [f64; 2])], rects: &Rects) -> Option> { + let (cu, cv) = (gx.len().saturating_sub(1), gy.len().saturating_sub(1)); + let cells = cu * cv; + let cell = |i: usize, j: usize| -> Option<([f64; 2], [f64; 2])> { + let (x0, x1, y0, y1) = (*gx.get(i)?, *gx.get(i + 1)?, *gy.get(j)?, *gy.get(j + 1)?); + // The clearance of the rectangle the cell belongs to, not of the (possibly thin) cell. + let b = rects.boxes.get(*rects.of.get(j * cu + i)?)?; + let (w, h) = (gx.get(b[2])? - gx.get(b[0])?, gy.get(b[3])? - gy.get(b[1])?); + let pad = CLEARANCE * w.min(h); + Some(([x0 - pad, y0 - pad], [x1 + pad, y1 + pad])) + }; + // Cells the boundary comes near: each segment marks the cells whose padded box it touches. + let mut near = vec![false; cells]; + let maxpad = rects + .boxes + .iter() + .filter_map(|b| Some(CLEARANCE * (gx.get(b[2])? - gx.get(b[0])?).min(gy.get(b[3])? - gy.get(b[1])?))) + .fold(0.0f64, f64::max); + let span = |g: &[f64], lo: f64, hi: f64| { + let a = g.partition_point(|t| *t < lo - maxpad).saturating_sub(1); + let b = (g.partition_point(|t| *t <= hi + maxpad) + 1).min(g.len().saturating_sub(1)); + (a, b) + }; + for (a, b) in segs { + let (i0, i1) = span(gx, a[0].min(b[0]), a[0].max(b[0])); + let (j0, j1) = span(gy, a[1].min(b[1]), a[1].max(b[1])); + for j in j0..j1 { + for i in i0..i1 { + let Some((lo, hi)) = cell(i, j) else { continue }; + if segment_meets_box(*a, *b, lo, hi) + && let Some(c) = near.get_mut(j * cu + i) + { + *c = true; + } + } + } + } + // The other cells are wholly inside or wholly outside: even–odd along each row's middle. + let mut keep = vec![true; rects.boxes.len()]; + for j in 0..cv { + let (Some(y0), Some(y1)) = (gy.get(j), gy.get(j + 1)) else { return None }; + let ym = 0.5 * (y0 + y1); + let mut cross: Vec = + segs.iter().filter(|(a, b)| (a[1] > ym) != (b[1] > ym)).map(|(a, b)| a[0] + (ym - a[1]) / (b[1] - a[1]) * (b[0] - a[0])).collect(); + cross.sort_by(f64::total_cmp); + for i in 0..cu { + let k = j * cu + i; + let (Some(x0), Some(x1)) = (gx.get(i), gx.get(i + 1)) else { return None }; + let xm = 0.5 * (x0 + x1); + let inside = cross.partition_point(|c| *c < xm) % 2 == 1; + if near.get(k).copied().unwrap_or(true) || !inside { + *keep.get_mut(*rects.of.get(k)?)? = false; + } + } + } + Some(keep) +} + +/// Whether the segment `a`–`b` meets the closed box `lo`–`hi` (Liang–Barsky clipping). +fn segment_meets_box(a: [f64; 2], b: [f64; 2], lo: [f64; 2], hi: [f64; 2]) -> bool { + let (mut t0, mut t1) = (0.0f64, 1.0f64); + for d in 0..2 { + let delta = b[d] - a[d]; + for (p, q) in [(-delta, a[d] - lo[d]), (delta, hi[d] - a[d])] { + if p == 0.0 { + if q < 0.0 { + return false; + } + } else { + let r = q / p; + if p < 0.0 { + t0 = t0.max(r); + } else { + t1 = t1.min(r); + } + if t0 > t1 { + return false; + } + } + } + } + true +} + +#[cfg(test)] +mod tests { + use super::*; + + fn circle(c: [f64; 2], r: f64, n: usize, ccw: bool) -> Vec<[f64; 2]> { + (0..n) + .map(|k| { + let a = std::f64::consts::TAU * k as f64 / n as f64 * if ccw { 1.0 } else { -1.0 }; + [c[0] + r * a.cos(), c[1] + r * a.sin()] + }) + .collect() + } + + /// A single loop over all `n` boundary points. + fn one(n: usize) -> [std::ops::Range; 1] { + [std::ops::Range { start: 0, end: n }] + } + + fn lines(a: f64, b: f64, n: usize) -> Vec { + (0..=n).map(|k| a + (b - a) * k as f64 / n as f64).collect() + } + + /// Every boundary segment is the side of exactly one triangle, every inner side of two, + /// and the triangles cover the region's area. + fn check(boundary: &[[f64; 2]], loops: &[std::ops::Range], area: f64, out: &(Vec<[f64; 2]>, Vec<[u32; 3]>)) { + let (extra, tris) = out; + let mut all = boundary.to_vec(); + all.extend_from_slice(extra); + let mut sides: std::collections::HashMap<(u32, u32), i32> = std::collections::HashMap::new(); + let mut sum = 0.0; + for t in tris { + let [a, b, c] = t.map(|i| all[i as usize]); + let ar = 0.5 * ((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])); + assert!(ar > 0.0, "counter-clockwise"); + sum += ar; + for e in 0..3 { + *sides.entry((t[e], t[(e + 1) % 3])).or_default() += 1; + } + } + for r in loops { + let n = r.len(); + for k in 0..n { + let (a, b) = ((r.start + k) as u32, (r.start + (k + 1) % n) as u32); + assert_eq!(sides.get(&(a, b)).copied().unwrap_or(0), 1, "boundary side {a}-{b} used once, forwards"); + sides.remove(&(a, b)); + } + } + for (&(a, b), &k) in &sides { + assert_eq!(k, 1); + assert_eq!(sides.get(&(b, a)).copied(), Some(1), "inner side {a}-{b} shared"); + } + assert!((sum - area).abs() < 1e-9 * area.max(1.0), "area {sum} vs {area}"); + } + + #[test] + fn a_disc_is_cells_and_a_band() { + let b = circle([0.0, 0.0], 1.0, 64, true); + let loops = one(64); + let poly_area = 0.5 * 64.0 * (std::f64::consts::TAU / 64.0).sin(); + let Some(out) = cut_grid(&b, &loops, &lines(-1.0, 1.0, 10), &lines(-1.0, 1.0, 10), &[], [1.0, 1.0]) else { panic!("meshes") }; + check(&b, &loops, poly_area, &out); + // Most of the disc is grid cells (two triangles each). + assert!(out.1.len() > 100, "{}", out.1.len()); + } + + #[test] + fn a_square_with_a_hole() { + let mut b = vec![[0.0, 0.0], [4.0, 0.0], [4.0, 4.0], [0.0, 4.0]]; + b.extend(circle([2.0, 2.0], 1.0, 32, false)); + let loops = [0..4, 4..36]; + let hole = 0.5 * 32.0 * (std::f64::consts::TAU / 32.0).sin(); + let Some(out) = cut_grid(&b, &loops, &lines(0.0, 4.0, 16), &lines(0.0, 4.0, 16), &[], [1.0, 1.0]) else { panic!("meshes") }; + check(&b, &loops, 16.0 - hole, &out); + } + + #[test] + fn a_curved_triangle_in_a_stretched_metric() { + // A three-cornered patch (as a vertex blend in its parameters), u twice as long as v. + let mut b: Vec<[f64; 2]> = (0..=8).map(|k| [k as f64 / 8.0, 0.0]).collect(); + b.extend((1..8).map(|k| { + let a = std::f64::consts::FRAC_PI_2 * k as f64 / 8.0; + [a.cos(), a.sin()] + })); + b.extend((1..8).rev().map(|k| [0.0, k as f64 / 8.0])); + let n = b.len(); + let loops = one(n); + let area: f64 = (0..n) + .map(|k| { + let (p, q) = (b[k], b[(k + 1) % n]); + 0.5 * (p[0] * q[1] - q[0] * p[1]) + }) + .sum(); + let Some(out) = cut_grid(&b, &loops, &lines(0.0, 1.0, 12), &lines(0.0, 1.0, 6), &[], [2.0, 1.0]) else { panic!("meshes") }; + check(&b, &loops, area, &out); + } + + /// The hole's motorcycle traces (lines through its four extremes) become T-mesh lines that + /// stop where the traces stop; the mesh stays exact on the boundary and covers the region. + #[test] + fn a_square_with_a_hole_and_its_traces() { + let mut b = vec![[0.0, 0.0], [4.0, 0.0], [4.0, 4.0], [0.0, 4.0]]; + b.extend(circle([1.7, 2.3], 1.0, 32, false)); + let loops = [0..4, 4..36]; + let hole = 0.5 * 32.0 * (std::f64::consts::TAU / 32.0).sin(); + let traces = super::super::motorcycle::traces(&b, &loops, [1.0, 1.0], 0.0); + assert_eq!(traces.len(), 4, "{traces:?}"); + let Some(out) = cut_grid(&b, &loops, &lines(0.0, 4.0, 16), &lines(0.0, 4.0, 16), &traces, [1.0, 1.0]) else { panic!("meshes") }; + check(&b, &loops, 16.0 - hole, &out); + // The traces run through the hole's extremes: along v = 2.3 from its u-extremes and + // along u = 1.7 from its v-extremes; their lines are in the grid. + for (d, x) in [(0, 1.7), (1, 2.3)] { + assert!(out.0.iter().any(|p| (p[d] - x).abs() < 1e-9), "a grid point on line {d} = {x}"); + } + } + + #[test] + fn nothing_kept_is_refused() { + let b = circle([0.0, 0.0], 0.1, 16, true); + assert!(cut_grid(&b, &one(16), &lines(-1.0, 1.0, 1), &lines(-1.0, 1.0, 1), &[], [1.0, 1.0]).is_none()); + } + + #[test] + fn segments_meet_boxes() { + assert!(segment_meets_box([-1.0, 0.5], [2.0, 0.5], [0.0, 0.0], [1.0, 1.0])); + assert!(!segment_meets_box([-1.0, 2.0], [2.0, 2.0], [0.0, 0.0], [1.0, 1.0])); + assert!(segment_meets_box([0.5, 0.5], [0.6, 0.6], [0.0, 0.0], [1.0, 1.0])); + assert!(!segment_meets_box([2.0, 0.0], [3.0, 1.0], [0.0, 0.0], [1.0, 1.0])); + } +} diff --git a/crates/kernel/src/tess/mod.rs b/crates/kernel/src/tess/mod.rs index 57ddce78..a5421b54 100644 --- a/crates/kernel/src/tess/mod.rs +++ b/crates/kernel/src/tess/mod.rs @@ -17,9 +17,13 @@ //! old way. pub(crate) mod cdt; +pub(crate) mod collar; pub(crate) mod coons; +pub(crate) mod cutgrid; pub(crate) mod edges; pub(crate) mod grid; +pub(crate) mod motorcycle; +pub(crate) mod pieces; pub(crate) mod size; pub(crate) mod ssi; pub(crate) mod surface; @@ -348,7 +352,7 @@ fn plan_face( let structured = match structured { None if ranges.len() == 1 => { let fine = |a: [f64; 2], b: [f64; 2]| uv::fine_between(s, rules, h, a, b); - coons::coons(&uvs, &fixed, speed, grid::MAX_FACE_GRID, &fine) + coons::coons(&uvs, &fixed, speed, grid::MAX_FACE_GRID, h, Some(&fine)) } other => other, }; @@ -362,12 +366,33 @@ fn plan_face( r } None => { - if planar { - stats.planar += 1; + // A curved face neither box nor four-cornered: cut into four-sided patches, else + // the iso-line grid cut along its motorcycle traces and trimmed by its boundary, + // else Delaunay. A plane: the butterfly of a round face, else Delaunay (a plane + // is shaded the same by any triangles, and its Delaunay needs no refinement). + let face = Region { s, uvs: &uvs, fixed: &fixed, ranges: &ranges, rules, h, speed, us: &us, vs: &vs }; + let trimmed = if planar { + if ranges.len() == 1 { plan_pieces(&face) } else { None } } else { - stats.delaunay += 1; + plan_pieces(&face).or_else(|| { + let traces = motorcycle::traces(&uvs, &ranges, speed, 0.25 * h); + cutgrid::cut_grid(&uvs, &ranges, &us, &vs, &traces, speed) + }) + }; + match trimmed { + Some(r) => { + stats.structured += 1; + r + } + None => { + if planar { + stats.planar += 1; + } else { + stats.delaunay += 1; + } + delaunay()? + } } - delaunay()? } } }; @@ -384,6 +409,211 @@ fn plan_face( Ok(FacePlan { face: u32::try_from(fi).unwrap_or(u32::MAX), surface: None, uv: all, fixed, normals: Vec::new(), flip: false, triangles }) } +/// A face's (u, v) region and what meshing it needs. +struct Region<'a> { + s: &'a EvalSurface, + uvs: &'a [[f64; 2]], + fixed: &'a [Vec3], + ranges: &'a [std::ops::Range], + rules: &'a size::Rules, + h: f64, + speed: [f64; 2], + us: &'a [f64], + vs: &'a [f64], +} + +/// A trimmed face cut into patches (squares round its holes or a butterfly in a round face, +/// joined to the outer loop by motorcycle traces), each meshed on its own (a Coons blend on its +/// best four corners; else the trimmed grid; else Delaunay), all sharing the cuts' nodes. The +/// cut nodes and patch interiors follow the boundary nodes. `None` when the cuts do not make +/// patches. +fn plan_pieces(f: &Region) -> Option<(Vec<[f64; 2]>, Vec<[u32; 3]>)> { + // The face's global seed: one spacing in every direction, the finest its size rules ask for + // along either parameter (the median gap between the grid lines they give, in mm). On a + // cylinder that is the spacing round it, which the cuts along its axis then get too. + let median_gap = |lines: &[f64], speed: f64| -> f64 { + let mut gaps: Vec = lines.windows(2).filter_map(|w| Some((w.get(1)? - w.first()?) * speed)).filter(|g| *g > 0.0).collect(); + gaps.sort_by(f64::total_cmp); + gaps.get(gaps.len() / 2).copied().unwrap_or(f64::INFINITY) + }; + let seed = median_gap(f.us, f.speed[0]).min(median_gap(f.vs, f.speed[1])).min(f.h); + let segments = collar::segments(f.uvs, f.ranges, f.speed, seed, matches!(f.s, EvalSurface::Plane { .. }))?; + let s = f.s; + let fine = |a: [f64; 2], b: [f64; 2]| uv::fine_between(s, f.rules, f.h, a, b); + // A cut's seeds: spaced from the boundary's own spacing at each end, growing geometrically + // between them, never wider than the global seed, and split further wherever the surface + // between two seeds fails the size rules. + let seeds = |a: [f64; 2], b: [f64; 2], sa: f64, sb: f64| -> Vec { graded(a, b, sa, sb, seed, f.speed, &fine) }; + let (cut_uv, pieces) = pieces::cut(f.uvs, f.ranges, &segments, f.speed, &seeds)?; + let n = f.uvs.len(); + let mut uv_all: Vec<[f64; 2]> = f.uvs.to_vec(); + uv_all.extend_from_slice(&cut_uv); + let mut pos: Vec = f.fixed.to_vec(); + for q in &cut_uv { + pos.push(s.derivs(q[0], q[1])?.0); + } + let mut extra: Vec<[f64; 2]> = cut_uv; + let mut tris: Vec<[u32; 3]> = Vec::new(); + for pc in &pieces { + // The patch's own loops, numbered from zero. + let mut ids: Vec = pc.ring.clone(); + let mut ranges = vec![std::ops::Range { start: 0, end: pc.ring.len() }]; + for hole in &pc.holes { + ranges.push(ids.len()..ids.len() + hole.len()); + ids.extend_from_slice(hole); + } + let luv: Vec<[f64; 2]> = ids.iter().map(|&k| uv_all.get(k as usize).copied()).collect::>()?; + let lpos: Vec = ids.iter().map(|&k| pos.get(k as usize).copied()).collect::>()?; + let mut bbox = [(f64::MAX, f64::MIN); 2]; + for p in &luv { + for d in 0..2 { + bbox[d] = (bbox[d].0.min(p[d]), bbox[d].1.max(p[d])); + } + } + let clip = |lines: &[f64], (a, b): (f64, f64)| -> Vec { + std::iter::once(a).chain(lines.iter().copied().filter(|x| *x > a && *x < b)).chain(std::iter::once(b)).collect() + }; + // A patch whose outline changes between boundary and cut at exactly four nodes (a + // quarter of a ring, a butterfly wing round a hole) takes those as its corners, so each + // side is all boundary or all cut; any other patch takes the best four of its corners. + let is_cut = |a: u32, b: u32| -> bool { + let (a, b) = (a as usize, b as usize); + a >= n || b >= n || !f.ranges.iter().any(|r| r.contains(&a) && r.contains(&b) && (a.abs_diff(b) == 1 || a.abs_diff(b) == r.len() - 1)) + }; + let ring = &pc.ring; + let switches: Vec = (0..ring.len()) + .filter(|&k| { + let (p, c, x) = (ring.get((k + ring.len() - 1) % ring.len()), ring.get(k), ring.get((k + 1) % ring.len())); + matches!((p, c, x), (Some(&p), Some(&c), Some(&x)) if is_cut(p, c) != is_cut(c, x)) + }) + .collect(); + // The patch is a rectangle in (u, v) with corners cut by free boundary (Lyon et al. §4; + // completed virtually in §4.2): its bounding box. A grid needs four real sides, so each + // corner is the node furthest out towards the box's corner (in the box's own + // proportions): an uncut corner is the corner itself, and on an arc cutting a corner it + // is the arc's 45° point, so the arc becomes a rounded corner shared by two sides. + let completed = || -> Option<[usize; 4]> { + let m = ring.len(); + let at = |k: usize| luv.get(k).map(|q| [q[0] * f.speed[0], q[1] * f.speed[1]]); + let (mut lo, mut hi) = ([f64::MAX; 2], [f64::MIN; 2]); + for k in 0..m { + let q = at(k)?; + for d in 0..2 { + lo[d] = lo[d].min(q[d]); + hi[d] = hi[d].max(q[d]); + } + } + let (w, h) = (hi[0] - lo[0], hi[1] - lo[1]); + if !(w > 0.0 && h > 0.0) { + return None; + } + let mut pick = [0usize; 4]; + for (slot, dir) in pick.iter_mut().zip([[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]]) { + *slot = (0..m).max_by(|&x, &y| { + let reach = |k: usize| at(k).map_or(f64::MIN, |q| dir[0] * (q[0] - lo[0]) / w + dir[1] * (q[1] - lo[1]) / h); + reach(x).total_cmp(&reach(y)) + })?; + } + // Four distinct corners in the loop's (counter-clockwise) order. + let start = pick[0]; + let rel = pick.map(|k| (k + m - start) % m); + (rel[0] < rel[1] && rel[1] < rel[2] && rel[2] < rel[3]).then_some(pick) + }; + let quad = if pc.holes.is_empty() { + // Every patch is four-sided by construction: its corners are where its outline turns. + coons::coons(&luv, &lpos, f.speed, grid::MAX_FACE_GRID, f.h, None).or_else(|| { + match switches[..] { + [a, b, c, d] => coons::coons_at(&luv, &lpos, [a, b, c, d], grid::MAX_FACE_GRID, f.h, None), + _ => None, + } + .or_else(|| completed().and_then(|c| coons::coons_at(&luv, &lpos, c, grid::MAX_FACE_GRID, f.h, None))) + }) + } else { + None + }; + let local = quad + .or_else(|| if pc.holes.is_empty() { coons::coons_best(&luv, &lpos, f.speed, grid::MAX_FACE_GRID, f.h) } else { None }) + .or_else(|| cutgrid::cut_grid(&luv, &ranges, &clip(f.us, bbox[0]), &clip(f.vs, bbox[1]), &[], f.speed)) + .or_else(|| { + let polys: Vec<&[[f64; 2]]> = ranges.iter().filter_map(|r| luv.get(r.clone())).collect(); + let seeds = cdt::lattice(bbox, [f.h / f.speed[0], f.h / f.speed[1]], &polys, f.speed); + let eval = |u: f64, v: f64| s.derivs(u, v).map(|(p, su, sv)| (p, su.cross(sv).normalized())); + let sag = cdt::Sag { eval: &eval, tol: f.rules.tol, h: f.h, turn: f.rules.turn, scale: f.speed }; + cdt::triangulate(&luv, &ranges, cdt::Fill { seeds: &seeds, refine: Some((sag, MAX_REFINE)), shape: true }) + })?; + let (lextra, ltris) = local; + let base = u32::try_from(n + extra.len()).ok()?; + let m = ids.len(); + for t in ltris { + let map = |i: u32| -> Option { + let i = i as usize; + if i < m { ids.get(i).copied() } else { base.checked_add(u32::try_from(i - m).ok()?) } + }; + tris.push([map(t[0])?, map(t[1])?, map(t[2])?]); + } + extra.extend(lextra); + } + Some((extra, tris)) +} + +/// Seeds on the segment `a`–`b` in (u, v) (fractions strictly between 0 and 1): node spacing +/// `sa` at `a` and `sb` at `b` (in the metric `speed`, each at most `h`), growing geometrically +/// between them (sizes `sa · q^k` summing to the length), then each interval split in half while +/// `fine` fails between its ends. +fn graded(a: [f64; 2], b: [f64; 2], sa: f64, sb: f64, h: f64, speed: [f64; 2], fine: &dyn Fn([f64; 2], [f64; 2]) -> bool) -> Vec { + let len = ((b[0] - a[0]) * speed[0]).hypot((b[1] - a[1]) * speed[1]); + let cap = |x: f64| { + if x.is_finite() && x > 0.0 { + x.min(if h > 0.0 { h } else { f64::INFINITY }) + } else if h > 0.0 { + h + } else { + len + } + }; + let (sa, sb) = (cap(sa), cap(sb)); + let mut cuts: Vec = Vec::new(); + if len.is_finite() && len > sa.max(sb) && sa > 0.0 && sb > 0.0 { + let q = (len - sa) / (len - sb); + let sizes: Vec = if (q - 1.0).abs() < 1e-3 || !(q > 0.0) || !q.is_finite() { + // Even spacing at the mean of the two. + let k = (len / (0.5 * (sa + sb))).round().clamp(1.0, 4096.0) as usize; + vec![len / k as f64; k] + } else { + let k = (1.0 + (sb / sa).ln() / q.ln()).round().clamp(1.0, 4096.0) as usize; + let raw: Vec = (0..k).map(|j| sa * q.powi(j as i32)).collect(); + let total: f64 = raw.iter().sum(); + raw.into_iter().map(|x| x * len / total).collect() + }; + let mut at = 0.0; + for s in sizes.iter().take(sizes.len().saturating_sub(1)) { + at += s; + cuts.push(at / len); + } + } + // The size rules between seeds (sag, facet angle). + let point = |t: f64| [a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t]; + let mut all: Vec = std::iter::once(0.0).chain(cuts).chain(std::iter::once(1.0)).collect(); + for _ in 0..12 { + let mut next = Vec::with_capacity(all.len() * 2); + let mut split = false; + for w in all.windows(2) { + let (Some(&x), Some(&y)) = (w.first(), w.get(1)) else { continue }; + next.push(x); + if !fine(point(x), point(y)) { + next.push(0.5 * (x + y)); + split = true; + } + } + next.push(1.0); + all = next; + if !split || all.len() > 4096 { + break; + } + } + all.into_iter().filter(|t| *t > 0.0 && *t < 1.0).collect() +} + /// Every boundary segment (consecutive nodes of a loop) is the side of exactly one triangle: /// the face covers its boundary with nothing missing or doubled, so it meets its neighbours /// exactly. Pole copies of one node (consecutive in a loop, at one point) are one vertex: the diff --git a/crates/kernel/src/tess/motorcycle.rs b/crates/kernel/src/tess/motorcycle.rs new file mode 100644 index 00000000..6dac8670 --- /dev/null +++ b/crates/kernel/src/tess/motorcycle.rs @@ -0,0 +1,354 @@ +//! The free-boundary motorcycle graph of a trimmed face (Lyon, Campen, Bommes & Kobbelt, +//! "Parametrization Quantization with Free Boundaries for Trimmed Quad Meshing", ACM TOG 38(4), +//! 2019, §4; after Eppstein, Goodrich, Kim & Tamstorf, "Motorcycle Graphs: Canonical Quad Mesh +//! Partitioning", SGP 2008). +//! +//! Peaks (Def. 4.1) are boundary nodes where u or v is a local extremum; a node in a flat run +//! counts only where the value differs on one side. A peak is concave where the face lies +//! beyond the extremum. Only concave peaks spawn motorcycles (Def. 4.2): a u-peak one along u, +//! a v-peak one along v, into the face. All move at once at the same speed (in the surface's +//! metric) along their iso-line; each stops where it reaches the boundary, a trace another +//! motorcycle laid earlier, or another motorcycle head-on (both stop where they meet). The +//! traces cut the face into patches that are rectangles with corners cut by free boundary. + +/// A trace: the iso-line `axis` = `at` (0: u = at, running along v; 1: v = at, running along u) +/// from `from` to `to` in the other parameter (`from < to`). +#[derive(Clone, Copy, Debug, PartialEq)] +pub(crate) struct Trace { + pub axis: usize, + pub at: f64, + pub from: f64, + pub to: f64, + /// The boundary node (peak) it starts from. + pub start: usize, + /// Its start and end points in (u, v). + pub a: [f64; 2], + pub b: [f64; 2], + /// Where it stopped. + pub end: End, +} + +/// Where a motorcycle stopped. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub(crate) enum End { + /// On the boundary segment between these two boundary nodes. + Boundary(usize, usize), + /// On the trace with this index (a T-node). + Trace(usize), + /// Head-on with the trace with this index (both end at one point). + HeadOn(usize), +} + +/// Most motorcycles one face spawns (the pairwise events grow with their square). +const MAX_MOTORCYCLES: usize = 512; + +struct Moto { + /// Origin (metric). + o: [f64; 2], + /// Direction of travel: the parameter (0: u, 1: v) and its sign. + d: usize, + sign: f64, + /// Boundary node it starts on (its own segments are not a crash). + node: usize, + /// Distance to the boundary along its way, and the segment it reaches there. + reach: f64, + hit: (usize, usize), +} + +/// The traces of the face bounded by `loops` (ranges of `boundary`, even–odd), in (u, v), +/// measured in the metric `scale`. Two motorcycles coming at each other on lines less than +/// `snap` apart (in the metric) meet head-on: a polygon's extreme node is only near the curve's +/// true extreme, so the two extremes facing each other across a gap (one hole above another) +/// are rarely on exactly one line. Empty when there is no concave peak. +pub(crate) fn traces(boundary: &[[f64; 2]], loops: &[std::ops::Range], scale: [f64; 2], snap: f64) -> Vec { + if !scale.iter().all(|s| s.is_finite() && *s > 0.0) { + return Vec::new(); + } + let p: Vec<[f64; 2]> = boundary.iter().map(|q| [q[0] * scale[0], q[1] * scale[1]]).collect(); + let mut lo = [f64::MAX; 2]; + let mut hi = [f64::MIN; 2]; + for q in &p { + for d in 0..2 { + lo[d] = lo[d].min(q[d]); + hi[d] = hi[d].max(q[d]); + } + } + let extent = (hi[0] - lo[0]).max(hi[1] - lo[1]); + if !(extent.is_finite() && extent > 0.0) { + return Vec::new(); + } + // Values this close are one (a node on an iso-line edge carries projection rounding). + let eps = extent * 1e-9; + let segs: Vec<(usize, usize)> = loops + .iter() + .flat_map(|r| { + let n = r.len(); + (0..n).map(move |k| (r.start + k, r.start + (k + 1) % n)) + }) + .collect(); + let inside = |q: [f64; 2]| { + let mut odd = false; + for &(a, b) in &segs { + let (Some(a), Some(b)) = (p.get(a), p.get(b)) else { continue }; + if (a[1] > q[1]) != (b[1] > q[1]) && q[0] < a[0] + (q[1] - a[1]) / (b[1] - a[1]) * (b[0] - a[0]) { + odd = !odd; + } + } + odd + }; + // Peaks and their motorcycles. + let mut motos: Vec = Vec::new(); + for r in loops { + let n = r.len(); + if n < 3 { + continue; + } + for k in 0..n { + let (Some(&a), Some(&b), Some(&c)) = (p.get(r.start + (k + n - 1) % n), p.get(r.start + k), p.get(r.start + (k + 1) % n)) else { + continue; + }; + for d in 0..2 { + let (before, after) = (a[d] - b[d], c[d] - b[d]); + let flat = |x: f64| x.abs() <= eps; + // A (non-strict) minimum higher on one side, or maximum lower on one side. + let min = (before > eps || flat(before)) && (after > eps || flat(after)) && !(flat(before) && flat(after)); + let max = (before < -eps || flat(before)) && (after < -eps || flat(after)) && !(flat(before) && flat(after)); + let sign = match (min, max) { + (true, false) => -1.0, + (false, true) => 1.0, + _ => continue, + }; + // Concave: the face lies beyond the extremum. + let mut probe = b; + probe[d] += sign * extent * 1e-6; + if !inside(probe) { + continue; + } + if motos.len() >= MAX_MOTORCYCLES { + return Vec::new(); + } + motos.push(Moto { o: b, d, sign, node: r.start + k, reach: f64::INFINITY, hit: (0, 0) }); + } + } + } + if motos.is_empty() { + return Vec::new(); + } + // Where each would reach the boundary. + for m in &mut motos { + let e = 1 - m.d; + for &(ia, ib) in &segs { + if ia == m.node || ib == m.node { + continue; + } + let (Some(a), Some(b)) = (p.get(ia), p.get(ib)) else { continue }; + // The segment crosses the line other = o[e] at parameter s along d. + let (lo_e, hi_e) = (a[e].min(b[e]), a[e].max(b[e])); + if m.o[e] < lo_e || m.o[e] > hi_e { + continue; + } + // A segment along the line itself is reached at its nearer end ahead. + let ahead = |x: f64| (x - m.o[m.d]) * m.sign; + let t = if (b[e] - a[e]).abs() <= eps { + [ahead(a[m.d]), ahead(b[m.d])].into_iter().filter(|t| *t > eps).fold(f64::INFINITY, f64::min) + } else { + ahead(a[m.d] + (m.o[e] - a[e]) / (b[e] - a[e]) * (b[m.d] - a[m.d])) + }; + if t > eps && t < m.reach { + m.reach = t; + m.hit = (ia, ib); + } + } + } + // Events in time order: (time, the motorcycle that stops, where, the one whose trace stops + // it and when it passed there). + struct Event { + t: f64, + who: usize, + at: f64, + cause: Option<(usize, f64)>, + partner: Option, + } + let mut events: Vec = Vec::new(); + for (i, m) in motos.iter().enumerate() { + if m.reach.is_finite() { + events.push(Event { t: m.reach, who: i, at: m.reach, cause: None, partner: None }); + } + for (j, n) in motos.iter().enumerate().skip(i + 1) { + if m.d != n.d { + // Crossing: m's line is other = m.o[1 - m.d]; n travels along n.d = 1 - m.d. + let x_on_m = (n.o[m.d] - m.o[m.d]) * m.sign; + let x_on_n = (m.o[n.d] - n.o[n.d]) * n.sign; + if x_on_m > eps && x_on_m < m.reach && x_on_n > eps && x_on_n < n.reach { + // The later one stops on the earlier one's trace (a tie stops both). + if x_on_m >= x_on_n { + events.push(Event { t: x_on_m, who: i, at: x_on_m, cause: Some((j, x_on_n)), partner: None }); + } + if x_on_n >= x_on_m { + events.push(Event { t: x_on_n, who: j, at: x_on_n, cause: Some((i, x_on_m)), partner: None }); + } + } + } else if m.sign != n.sign && (m.o[1 - m.d] - n.o[1 - n.d]).abs() <= eps.max(snap) { + // Head-on on one line: they meet halfway if each would get there. + let gap = (n.o[m.d] - m.o[m.d]) * m.sign; + let half = 0.5 * gap; + if gap > eps && half < m.reach && half < n.reach { + events.push(Event { t: half, who: i, at: half, cause: None, partner: Some(j) }); + } + } + } + } + events.sort_by(|a, b| a.t.total_cmp(&b.t)); + let mut stop: Vec> = vec![None; motos.len()]; + for ev in &events { + if stop.get(ev.who).is_none_or(|s| s.is_some()) { + continue; + } + if let Some((c, passed)) = ev.cause { + // The trace is there only if its motorcycle got that far. + if stop.get(c).copied().flatten().is_some_and(|(s, _)| s < passed) { + continue; + } + } + let how = match (ev.cause, ev.partner) { + (Some((c, _)), _) => Stop::Trace(c), + (None, Some(j)) => Stop::HeadOn(j), + (None, None) => Stop::Boundary, + }; + if let Some(j) = ev.partner { + if stop.get(j).is_none_or(|s| s.is_some()) { + continue; + } + if let Some(s) = stop.get_mut(j) { + *s = Some((ev.at, Stop::HeadOn(ev.who))); + } + } + if let Some(s) = stop.get_mut(ev.who) { + *s = Some((ev.at, how)); + } + } + // Output indices of the motorcycles that stopped (all do: each reaches the boundary). + let mut out_index = vec![usize::MAX; motos.len()]; + let mut k = 0; + for (i, s) in stop.iter().enumerate() { + if s.is_some() { + if let Some(slot) = out_index.get_mut(i) { + *slot = k; + } + k += 1; + } + } + let uv = |q: [f64; 2]| [q[0] / scale[0], q[1] / scale[1]]; + motos + .iter() + .zip(&stop) + .filter_map(|(m, s)| { + let (len, how) = (*s)?; + let (a, b) = (m.o[m.d], m.o[m.d] + m.sign * len); + let (from, to) = (a.min(b) / scale[m.d], a.max(b) / scale[m.d]); + // The trace runs along d, so it lies on the iso-line of the other parameter. + let axis = 1 - m.d; + let mut endp = m.o; + endp[m.d] = b; + let end = match how { + Stop::Boundary => End::Boundary(m.hit.0, m.hit.1), + Stop::Trace(c) => End::Trace(*out_index.get(c)?), + Stop::HeadOn(j) => End::HeadOn(*out_index.get(j)?), + }; + Some(Trace { axis, at: m.o[axis] / scale[axis], from, to, start: m.node, a: uv(m.o), b: uv(endp), end }) + }) + .collect() +} + +#[derive(Clone, Copy, Debug)] +enum Stop { + Boundary, + Trace(usize), + HeadOn(usize), +} + +#[cfg(test)] +mod tests { + use super::*; + + fn circle(c: [f64; 2], r: f64, n: usize, ccw: bool) -> Vec<[f64; 2]> { + (0..n) + .map(|k| { + let a = std::f64::consts::TAU * k as f64 / n as f64 * if ccw { 1.0 } else { -1.0 }; + [c[0] + r * a.cos(), c[1] + r * a.sin()] + }) + .collect() + } + + #[test] + fn a_convex_face_has_no_traces() { + let b = circle([0.0, 0.0], 1.0, 32, true); + assert!(traces(&b, &[std::ops::Range { start: 0, end: 32 }], [1.0, 1.0], 0.0).is_empty()); + let sq = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]; + assert!(traces(&sq, &[std::ops::Range { start: 0, end: 4 }], [1.0, 1.0], 0.0).is_empty()); + } + + /// A centred round hole: its four extremes are concave peaks, and each motorcycle runs + /// straight to the outer square (none meets another). + #[test] + fn a_hole_sends_four_traces_to_the_outside() { + let mut b = vec![[0.0, 0.0], [4.0, 0.0], [4.0, 4.0], [0.0, 4.0]]; + b.extend(circle([2.0, 2.0], 1.0, 32, false)); + let t = traces(&b, &[0..4, 4..36], [1.0, 1.0], 0.0); + assert_eq!(t.len(), 4, "{t:?}"); + let on = |axis: usize, at: f64, from: f64, to: f64| { + t.iter().any(|x| x.axis == axis && (x.at - at).abs() < 1e-12 && (x.from - from).abs() < 1e-9 && (x.to - to).abs() < 1e-9) + }; + assert!(on(1, 2.0, 0.0, 1.0), "left: v = 2 from u = 0 to the hole"); + assert!(on(1, 2.0, 3.0, 4.0), "right"); + assert!(on(0, 2.0, 0.0, 1.0), "below"); + assert!(on(0, 2.0, 3.0, 4.0), "above"); + } + + /// An L: the reflex corner is a concave peak in both u and v, and its two motorcycles + /// cut the L into three rectangles. + #[test] + fn an_l_is_cut_at_its_reflex_corner() { + let b = [[0.0, 0.0], [2.0, 0.0], [2.0, 1.0], [1.0, 1.0], [1.0, 2.0], [0.0, 2.0]]; + let t = traces(&b, &[std::ops::Range { start: 0, end: 6 }], [1.0, 1.0], 0.0); + assert_eq!(t.len(), 2, "{t:?}"); + let has = |axis: usize| { + t.iter().any(|x| x.axis == axis && x.at == 1.0 && x.from == 0.0 && x.to == 1.0 && x.start == 3 && matches!(x.end, End::Boundary(..))) + }; + assert!(has(0) && has(1), "{t:?}"); + } + + /// Two holes side by side: the inner motorcycles meet head-on halfway, and a later + /// motorcycle stops on an earlier one's trace (no crossings). + #[test] + fn traces_meet_head_on_and_never_cross() { + let mut b = vec![[0.0, 0.0], [8.0, 0.0], [8.0, 4.0], [0.0, 4.0]]; + b.extend(circle([2.0, 2.0], 1.0, 32, false)); + b.extend(circle([6.0, 2.5], 1.0, 32, false)); + let t = traces(&b, &[0..4, 4..36, 36..68], [1.0, 1.0], 0.0); + // Hole 1's right extreme (3, 2) and hole 2's left extreme (5, 2.5) are on different + // lines; each runs on until it meets the other hole or a trace. + for (i, a) in t.iter().enumerate() { + for c in t.iter().skip(i + 1) { + if a.axis != c.axis { + // a on axis line a.at; c crosses it if c.at is strictly within a's span + // and a.at strictly within c's. + let cross = c.at > a.from + 1e-9 && c.at < a.to - 1e-9 && a.at > c.from + 1e-9 && a.at < c.to - 1e-9; + assert!(!cross, "{a:?} crosses {c:?}"); + } + } + } + // Same-line pair: two holes at the same height meet halfway. + let mut b = vec![[0.0, 0.0], [8.0, 0.0], [8.0, 4.0], [0.0, 4.0]]; + b.extend(circle([2.0, 2.0], 1.0, 32, false)); + b.extend(circle([6.0, 2.0], 1.0, 32, false)); + let t = traces(&b, &[0..4, 4..36, 36..68], [1.0, 1.0], 0.0); + assert!(t.iter().any(|x| x.axis == 1 && (x.at - 2.0).abs() < 1e-12 && (x.from - 3.0).abs() < 1e-9 && (x.to - 4.0).abs() < 1e-9), "{t:?}"); + assert!(t.iter().any(|x| x.axis == 1 && (x.at - 2.0).abs() < 1e-12 && (x.from - 4.0).abs() < 1e-9 && (x.to - 5.0).abs() < 1e-9), "{t:?}"); + // Those two ended head-on with each other. + let i = t.iter().position(|x| x.axis == 1 && (x.from - 3.0).abs() < 1e-9).unwrap_or(usize::MAX); + let j = t.iter().position(|x| x.axis == 1 && (x.from - 4.0).abs() < 1e-9).unwrap_or(usize::MAX); + assert_eq!(t.get(i).map(|x| x.end), Some(End::HeadOn(j)), "{t:?}"); + assert_eq!(t.get(j).map(|x| x.end), Some(End::HeadOn(i)), "{t:?}"); + } +} diff --git a/crates/kernel/src/tess/pieces.rs b/crates/kernel/src/tess/pieces.rs new file mode 100644 index 00000000..9aa76273 --- /dev/null +++ b/crates/kernel/src/tess/pieces.rs @@ -0,0 +1,348 @@ +//! A trimmed face cut into four-sided patches. The cuts are segments in (u, v): the squares +//! round the face's holes (each joined to its hole at four points, so the ring between them is +//! four collar patches) and the motorcycle traces that extend the squares' sides. Each cut is +//! sampled once and both patches beside it use those nodes (as faces share their edges' nodes), +//! so the patches meet exactly. A cut never adds a node to the face's boundary: it ends on a +//! boundary node, at a point inside the face (a square's corner, or where two traces met +//! head-on), or on another cut (a T-node, which that cut then passes through). + +use std::collections::HashMap; + +/// Where a cut segment ends. +#[derive(Clone, Copy, Debug, PartialEq)] +pub(crate) enum End { + /// A boundary node. + Node(u32), + /// A point inside the face; points at the same place are one node. + Point([f64; 2]), + /// A point on cut segment `.0` (a T-node: that segment passes through it). + On(usize, [f64; 2]), +} + +/// A cut from `a` to `b`, straight in (u, v). +#[derive(Clone, Copy, Debug, PartialEq)] +pub(crate) struct Segment { + pub a: End, + pub b: End, + /// Segments to sample it in (spread over its pieces between T-nodes by length), instead of + /// the seed rule: a butterfly's square side takes the count of the hole arc facing it. + pub count: Option, +} + +/// A patch: its outer ring and any holes inside it, as indices into the boundary nodes followed +/// by the cut nodes, counter-clockwise (holes clockwise) in (u, v). +#[derive(Clone, Debug, PartialEq, Eq)] +pub(crate) struct Piece { + pub ring: Vec, + pub holes: Vec>, +} + +/// Cut the face bounded by `loops` (ranges of `boundary`) along `segments`, measured in the +/// metric `scale`. `seeds(a, b, sa, sb)`: where to put nodes between two nodes of a cut (in +/// (u, v)), as fractions strictly between 0 and 1, given the spacing wanted at each (in the +/// metric): the boundary's own spacing at a boundary node, unbounded inside the face (the seed +/// rule then sets it). Returns the cut nodes (numbered after the boundary) and the patches; +/// `None` when the cuts do not make a planar subdivision (the face is then meshed whole). +pub(crate) fn cut( + boundary: &[[f64; 2]], + loops: &[std::ops::Range], + segments: &[Segment], + scale: [f64; 2], + seeds: &dyn Fn([f64; 2], [f64; 2], f64, f64) -> Vec, +) -> Option<(Vec<[f64; 2]>, Vec)> { + let n = boundary.len(); + let metric = |q: [f64; 2]| [q[0] * scale[0], q[1] * scale[1]]; + let dist = |a: [f64; 2], b: [f64; 2]| { + let (a, b) = (metric(a), metric(b)); + (a[0] - b[0]).hypot(a[1] - b[1]) + }; + let mut cut_uv: Vec<[f64; 2]> = Vec::new(); + let mut points: HashMap<[u64; 2], u32> = HashMap::new(); + // Each segment's end nodes, and the T-nodes other segments put on it. + let mut ends: Vec<(u32, u32)> = Vec::with_capacity(segments.len()); + let mut along: Vec> = vec![Vec::new(); segments.len()]; + for s in segments { + let mut key = |e: End| -> Option { + let mut point = |q: [f64; 2]| -> Option { + let bits = [q[0].to_bits(), q[1].to_bits()]; + if let Some(&k) = points.get(&bits) { + return Some(k); + } + let k = u32::try_from(n + cut_uv.len()).ok()?; + cut_uv.push(q); + points.insert(bits, k); + Some(k) + }; + match e { + End::Node(k) => ((k as usize) < n).then_some(k), + End::Point(q) => point(q), + End::On(on, q) => { + let k = point(q)?; + along.get_mut(on)?.push(k); + Some(k) + } + } + }; + let (a, b) = (key(s.a)?, key(s.b)?); + ends.push((a, b)); + } + let at_in = |cut_uv: &[[f64; 2]], k: u32| -> Option<[f64; 2]> { + let k = k as usize; + if k < n { boundary.get(k).copied() } else { cut_uv.get(k - n).copied() } + }; + // The spacing wanted at a node: at a boundary node the mean of its two boundary segments. + let spacing = |k: u32| -> f64 { + let k = k as usize; + let Some(r) = loops.iter().find(|r| r.contains(&k)) else { return f64::INFINITY }; + let m = r.len(); + let i = k - r.start; + let (prev, next) = (r.start + (i + m - 1) % m, r.start + (i + 1) % m); + match (boundary.get(prev), boundary.get(k), boundary.get(next)) { + (Some(&p), Some(&q), Some(&x)) => 0.5 * (dist(p, q) + dist(q, x)), + _ => f64::INFINITY, + } + }; + // Boundary sides, so a cut that is one of them is left out. + let mut sides: std::collections::HashSet<(u32, u32)> = std::collections::HashSet::new(); + for r in loops { + let m = r.len(); + for k in 0..m { + let (a, b) = (u32::try_from(r.start + k).ok()?, u32::try_from(r.start + (k + 1) % m).ok()?); + sides.insert((a.min(b), a.max(b))); + } + } + let mut cuts: Vec<(u32, u32)> = Vec::new(); + let mut taken: std::collections::HashSet<(u32, u32)> = std::collections::HashSet::new(); + for (i, &(s0, s1)) in ends.iter().enumerate() { + if s0 == s1 { + continue; + } + let p0 = at_in(&cut_uv, s0)?; + let mut nodes: Vec<(f64, u32)> = vec![(0.0, s0)]; + for &k in along.get(i)? { + nodes.push((dist(p0, at_in(&cut_uv, k)?), k)); + } + let total = dist(p0, at_in(&cut_uv, s1)?); + nodes.push((total, s1)); + nodes.sort_by(|a, b| a.0.total_cmp(&b.0)); + nodes.dedup_by_key(|x| x.1); + let count = segments.get(i).and_then(|s| s.count); + for w in nodes.windows(2) { + let (Some(&(_, a)), Some(&(_, b))) = (w.first(), w.get(1)) else { continue }; + if a == b { + continue; + } + // The same cut twice (two segments between the same nodes) is one cut. + let key = (a.min(b), a.max(b)); + if sides.contains(&key) || !taken.insert(key) { + continue; + } + let (pa, pb) = (at_in(&cut_uv, a)?, at_in(&cut_uv, b)?); + let fractions: Vec = match count { + Some(k) if total > 0.0 => { + let share = dist(pa, pb) / total; + let k = ((k as f64 * share).round() as usize).clamp(1, 4096); + (1..k).map(|j| j as f64 / k as f64).collect() + } + _ => seeds(pa, pb, spacing(a), spacing(b)), + }; + let mut prev = a; + for f in fractions.into_iter().take(4096) { + let k = u32::try_from(n + cut_uv.len()).ok()?; + cut_uv.push([pa[0] + (pb[0] - pa[0]) * f, pa[1] + (pb[1] - pa[1]) * f]); + cuts.push((prev, k)); + prev = k; + } + cuts.push((prev, b)); + } + } + if cuts.is_empty() { + return None; + } + let at = |k: u32| at_in(&cut_uv, k); + let pos = |k: u32| at(k).map(metric); + // Half-edges: boundary sides with the face on their left, cut sides both ways. + let inside = |q: [f64; 2]| { + let mut odd = false; + for r in loops { + let m = r.len(); + for k in 0..m { + let (Some(a), Some(b)) = (boundary.get(r.start + k).map(|p| metric(*p)), boundary.get(r.start + (k + 1) % m).map(|p| metric(*p))) + else { + continue; + }; + if (a[1] > q[1]) != (b[1] > q[1]) && q[0] < a[0] + (q[1] - a[1]) / (b[1] - a[1]) * (b[0] - a[0]) { + odd = !odd; + } + } + } + odd + }; + let mut half: Vec<(u32, u32)> = Vec::new(); + for r in loops { + let m = r.len(); + let side = |k: usize| -> Option<(u32, u32)> { Some((u32::try_from(r.start + k).ok()?, u32::try_from(r.start + (k + 1) % m).ok()?)) }; + // Which way round the face lies to the left: test beside the loop's longest side. + let longest = (0..m).max_by(|&x, &y| { + let l = |k: usize| side(k).and_then(|(a, b)| Some(dist(at(a)?, at(b)?))).unwrap_or(0.0); + l(x).total_cmp(&l(y)) + })?; + let (a, b) = side(longest)?; + let (pa, pb) = (pos(a)?, pos(b)?); + let len = (pb[0] - pa[0]).hypot(pb[1] - pa[1]); + if !(len > 0.0) { + return None; + } + let left = [0.5 * (pa[0] + pb[0]) - (pb[1] - pa[1]) * 1e-4, 0.5 * (pa[1] + pb[1]) + (pb[0] - pa[0]) * 1e-4]; + let forward = inside(left); + for k in 0..m { + let (a, b) = side(k)?; + if a == b { + continue; + } + half.push(if forward { (a, b) } else { (b, a) }); + } + } + for &(a, b) in &cuts { + half.push((a, b)); + half.push((b, a)); + } + // Outgoing half-edges per node, by angle. + let angle = |a: u32, b: u32| -> Option { + let (p, q) = (pos(a)?, pos(b)?); + Some((q[1] - p[1]).atan2(q[0] - p[0])) + }; + let mut out: HashMap> = HashMap::new(); + for (h, &(a, b)) in half.iter().enumerate() { + out.entry(a).or_default().push((angle(a, b)?, h)); + } + for v in out.values_mut() { + v.sort_by(|x, y| x.0.total_cmp(&y.0)); + } + // The face on the left of u→v continues at v with the first outgoing side clockwise from + // v→u. + let next = |h: usize| -> Option { + let &(u, v) = half.get(h)?; + let back = angle(v, u)?; + let list = out.get(&v)?; + let before = list + .iter() + .rev() + .find(|(a, e)| *a < back - 1e-12 && half.get(*e).is_some_and(|x| x.1 != u)) + .or_else(|| list.iter().rev().find(|(_, e)| half.get(*e).is_some_and(|x| x.1 != u))); + before.map(|(_, e)| *e).or_else(|| list.iter().find(|(_, e)| half.get(*e).is_some_and(|x| x.1 == u)).map(|(_, e)| *e)) + }; + let mut seen = vec![false; half.len()]; + let mut rings: Vec> = Vec::new(); + for h0 in 0..half.len() { + if *seen.get(h0)? { + continue; + } + let mut ring = Vec::new(); + let mut h = h0; + loop { + *seen.get_mut(h)? = true; + ring.push(half.get(h)?.0); + h = next(h)?; + if h == h0 { + break; + } + if ring.len() > half.len() || *seen.get(h)? { + return None; + } + } + rings.push(ring); + } + let area = |ring: &[u32]| -> Option { + let mut s = 0.0; + for k in 0..ring.len() { + let (a, b) = (pos(*ring.get(k)?)?, pos(*ring.get((k + 1) % ring.len())?)?); + s += a[0] * b[1] - b[0] * a[1]; + } + Some(0.5 * s) + }; + let mut pieces: Vec = Vec::new(); + let mut holes: Vec> = Vec::new(); + for r in rings { + if area(&r)? > 0.0 { + pieces.push(Piece { ring: r, holes: Vec::new() }); + } else { + holes.push(r); + } + } + // A loop no cut reached is a hole of the patch around it. + for h in holes { + let p = pos(*h.first()?)?; + let contains = |ring: &[u32]| { + let mut odd = false; + for k in 0..ring.len() { + let (Some(a), Some(b)) = (ring.get(k).and_then(|x| pos(*x)), ring.get((k + 1) % ring.len()).and_then(|x| pos(*x))) else { continue }; + if (a[1] > p[1]) != (b[1] > p[1]) && p[0] < a[0] + (p[1] - a[1]) / (b[1] - a[1]) * (b[0] - a[0]) { + odd = !odd; + } + } + odd + }; + let host = pieces.iter_mut().filter(|pc| contains(&pc.ring)).min_by(|a, b| { + let s = |pc: &Piece| area(&pc.ring).unwrap_or(f64::MAX); + s(a).total_cmp(&s(b)) + })?; + host.holes.push(h); + } + Some((cut_uv, pieces)) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn area(all: &[[f64; 2]], ring: &[u32]) -> f64 { + (0..ring.len()) + .map(|k| { + let (a, b) = (all[ring[k] as usize], all[ring[(k + 1) % ring.len()] as usize]); + 0.5 * (a[0] * b[1] - b[0] * a[1]) + }) + .sum() + } + + /// A square cut by four segments meeting at its centre (one point inside the face, shared) + /// gives four quarters covering it. + #[test] + fn a_cross_cuts_a_square_into_four() { + let b = vec![[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [0.0, 2.0], [0.0, 1.0]]; + let loops = [std::ops::Range { start: 0, end: 8 }]; + let c = End::Point([1.0, 1.0]); + let segs = [ + Segment { a: End::Node(1), b: c, count: None }, + Segment { a: End::Node(3), b: c, count: None }, + Segment { a: End::Node(5), b: c, count: None }, + Segment { a: End::Node(7), b: c, count: None }, + ]; + let Some((cut_uv, pieces)) = cut(&b, &loops, &segs, [1.0, 1.0], &|_, _, _, _| vec![0.5]) else { panic!("cuts") }; + assert_eq!(pieces.len(), 4); + let mut all = b.clone(); + all.extend(cut_uv); + let total: f64 = pieces.iter().map(|p| area(&all, &p.ring)).sum(); + assert!((total - 4.0).abs() < 1e-12, "{total}"); + } + + /// A segment ending on another (a T-node) splits it there: three patches. + #[test] + fn a_t_node_splits_the_segment_it_lands_on() { + let b = vec![[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [2.0, 1.0], [2.0, 2.0], [0.0, 2.0], [0.0, 1.0]]; + let loops = [std::ops::Range { start: 0, end: 7 }]; + let segs = [Segment { a: End::Node(6), b: End::Node(3), count: None }, Segment { a: End::Node(1), b: End::On(0, [1.0, 1.0]), count: None }]; + let Some((cut_uv, pieces)) = cut(&b, &loops, &segs, [1.0, 1.0], &|_, _, _, _| Vec::new()) else { panic!("cuts") }; + assert_eq!(pieces.len(), 3, "{pieces:?}"); + let mut all = b.clone(); + all.extend(cut_uv); + let total: f64 = pieces.iter().map(|p| area(&all, &p.ring)).sum(); + assert!((total - 4.0).abs() < 1e-12, "{total}"); + } + + #[test] + fn no_segments_no_cut() { + let b = [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]; + assert!(cut(&b, &[std::ops::Range { start: 0, end: 4 }], &[], [1.0, 1.0], &|_, _, _, _| Vec::new()).is_none()); + } +} diff --git a/crates/kernel/src/tess/size.rs b/crates/kernel/src/tess/size.rs index 37129878..9cbad54c 100644 --- a/crates/kernel/src/tess/size.rs +++ b/crates/kernel/src/tess/size.rs @@ -14,10 +14,9 @@ use super::uv::Domain; /// Longest side, as a fraction of the body's extent. const MAX_LEN_FRACTION: f64 = 1.0 / 5.0; /// Longest side relative to the side the curvature asks for: along a straight direction (a -/// cylinder's axis) triangles may be this many times longer than across it. With a 30° facet -/// angle and [`MAX_LEN_FRACTION`] this puts the oracle parts' display meshes at about the -/// reference mesher's triangle count (1.2×; 5× with equilateral triangles and 15°). -const ASPECT: f64 = 4.0; +/// cylinder's axis) cells may be this many times longer than across it. 2: cells stay close to +/// square (patches mesh as grids) without doubling the rows along every straight surface. +const ASPECT: f64 = 2.0; /// Shortest side the aspect limit may ask for, as a fraction of the longest side. A bend whose /// radius is tiny next to the body (a hem) would otherwise tie the whole length of its /// straight edges to that radius; along a straight direction the normal does not turn, so @@ -172,15 +171,16 @@ mod tests { transform: [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]], swap: false, }; - // A 20 mm body: sides from 0.5 mm (the floor) to 4 mm. + // A 20 mm body: sides from 1 mm (the floor, a quarter of the longest) to 4 mm. let rules = Rules::new(0.01, 20.0, THIRTY); let (h, speed) = face_size(&s, &super::super::uv::domain(&s), &rules); - // r = 3: 30° of turn is 1.57 mm; the 0.01 mm sag chord is √(8·3·0.01) = 0.49 mm. - assert!((h - ASPECT * (8.0 * 3.0 * 0.01f64).sqrt()).abs() < 0.01, "{h}"); + // r = 3: 30° of turn is 1.57 mm; the 0.01 mm sag chord is √(8·3·0.01) = 0.49 mm, times the + // aspect, unless that is under the floor. + assert!((h - (ASPECT * (8.0 * 3.0 * 0.01f64).sqrt()).max(20.0 * MAX_LEN_FRACTION * MIN_LEN_FRACTION)).abs() < 0.01, "{h}"); assert!((speed[0] - 40.0).abs() < 1e-6 && (speed[1] - 3.0).abs() < 1e-6, "{speed:?}"); // Coarser: the turn limit governs. let (h, _) = face_size(&s, &super::super::uv::domain(&s), &Rules::new(0.5, 100.0, THIRTY)); - assert!((h - ASPECT * 3.0 * std::f64::consts::PI / 6.0).abs() < 0.01, "{h}"); + assert!((h - (ASPECT * 3.0 * std::f64::consts::PI / 6.0).max(100.0 * MAX_LEN_FRACTION * MIN_LEN_FRACTION)).abs() < 0.01, "{h}"); // Without a facet angle the sag alone: √(8·3·0.5) = 3.46 mm, four times along. let (h, _) = face_size(&s, &super::super::uv::domain(&s), &Rules::new(0.5, 100.0, None)); assert!((h - ASPECT * (8.0 * 3.0 * 0.5f64).sqrt()).abs() < 0.01, "{h}"); diff --git a/crates/kernel/src/tess/tests.rs b/crates/kernel/src/tess/tests.rs index 6580299c..5587fb7d 100644 --- a/crates/kernel/src/tess/tests.rs +++ b/crates/kernel/src/tess/tests.rs @@ -108,13 +108,14 @@ fn holes_and_trimmed_faces_are_watertight() { let bx = ok(crate::box_solid(Vec3::ZERO, Vec3::new(10.0, 10.0, 4.0))); let pin = ok(crate::cylinder(Vec3::new(5.0, 5.0, -1.0), Vec3::Z, 2.0, 6.0)); check("box with a hole", &cut(&bx, &pin), tol, 400.0 - PI * 4.0 * 4.0); - // A cylinder with a cross hole: trimmed cylinder halves (Delaunay) and intersection-curve walls. + // A cylinder with a cross hole: trimmed cylinder halves (an iso-line grid cut by the hole) + // and intersection-curve walls; no face needs Delaunay. let cyl = ok(crate::cylinder(Vec3::new(5.0, 5.0, -1.0), Vec3::Z, 2.0, 6.0)); let cross = ok(crate::cylinder(Vec3::new(5.0, -1.0, 2.0), Vec3::Y, 1.0, 12.0)); let body = cut(&cyl, &cross); let (m, stats) = plan_body(&body, tol); assert_watertight(&m, "cross hole"); - assert!(stats.delaunay > 0, "{stats:?}"); + assert_eq!(stats.delaunay, 0, "{stats:?}"); // Against the reference mesher's volume (no closed form for the two-cylinder intersection). let Ok(reference) = body.tessellate(tol * 0.2) else { panic!("reference mesh") }; let (a, b) = (m.measure().volume, reference.measure().volume); @@ -279,3 +280,19 @@ fn a_placed_body_meshes_like_the_original() { let (p, _) = plan_body(&placed, 0.005); assert_eq!(m.tri_face, p.tri_face); } + +/// Seeds on a motorcycle trace start at the spacing of the boundary it leaves and grow +/// geometrically to the spacing at its other end (here 0.1 to 1 over 5 mm). +#[test] +fn trace_seeds_grade_from_one_end_to_the_other() { + let t = graded([0.0, 0.0], [5.0, 0.0], 0.1, 1.0, 10.0, [1.0, 1.0], &|_, _| true); + let all: Vec = std::iter::once(0.0).chain(t.iter().copied()).chain(std::iter::once(1.0)).map(|f| f * 5.0).collect(); + let sizes: Vec = all.windows(2).map(|w| w[1] - w[0]).collect(); + assert!((sizes[0] - 0.1).abs() < 0.03, "first {}", sizes[0]); + assert!((sizes[sizes.len() - 1] - 1.0).abs() < 0.3, "last {}", sizes[sizes.len() - 1]); + assert!(sizes.windows(2).all(|w| w[1] >= w[0] - 1e-12), "growing: {sizes:?}"); + // Equal ends: even; a segment no longer than its spacing: no seeds. + let even = graded([0.0, 0.0], [4.0, 0.0], 1.0, 1.0, 10.0, [1.0, 1.0], &|_, _| true); + assert_eq!(even.len(), 3); + assert!(graded([0.0, 0.0], [0.5, 0.0], 1.0, 1.0, 10.0, [1.0, 1.0], &|_, _| true).is_empty()); +}