From 201e9d5f9f0332da9b9bdec31d175dec71bbaf70 Mon Sep 17 00:00:00 2001 From: Irakli Kourtis Date: Mon, 12 Oct 2026 21:37:01 -0400 Subject: [PATCH] Patch fills non-planar edge loops with one smooth face (#250) (#472) * Add a one-surface fill of closed 3D loops to geom A closed loop of curves that does not lie in a plane has no exact surface; the usual answer is one smooth face fitted over it. fill_loop does that: - the loop's best-fit plane (Newell normal), with a frame chosen from the loop itself so the result does not depend on where the loop starts, its direction, or where the part sits - the loop seen along that normal must be simple and must not fold or run along the normal; otherwise it is refused ("not supported yet"), never filled with wrong geometry - one non-rational bicubic B-spline height field over the loop's box, fitted by least squares to the loop's samples (The NURBS Book 9.4.1, 9.4.3, eq. 9.62/9.65) with a thin-plate smoothness term lambda*A*int(h_uu^2 + 2 h_uv^2 + h_vv^2), A one grid cell's area, so the fit is the same at any scale; only planes are reproduced exactly - span counts per direction (long narrow loops fit), refined until the loop is within tolerance, else refused with the deviation reached bspline: bisected span search, basis functions on a stack array, one clamping helper shared by curves and surfaces, one binomial helper. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01J53sR6gEaVt56soFDpTck8 * Patch fills non-planar loops of edges with one smooth face (#250) surface.patch with edges refused any loop that is not planar ("not supported yet: patching a loop that is not planar"). It now fills it: - planar loops still get one exact planar face, but planarity is certified conservatively: lines by their ends, B-spline/NURBS edges by their control points, intersection curves by their samples lying on a support plane within tolerance; a loop that can't be certified is filled instead - other loops get one face from geom's fill_loop, bounded by the body's own chained edges, whole (no split): its history is one generated face whose edges are the body's, and it stitches back into a closed solid - intersection-curve edges are sampled by a seeded, guarded double projection shared by the planarity test and the fill (plane_fit) - the face tolerance is 10x the planarity tolerance (documented: later operations take the face to contain its edges within it) - edge picking runs inside the guard; the geom->truck conversion is shared with the T-spline sheets (spline_conv); the Newell plane fit is shared with STEP import Tests: 3- to 8-sided, concave (dart, arrow, plus), flat-corner, two-edge and imported rims fill as one face and stitch closed (volumes against an independent integral); finite mesh normals, same result for any pick order or placement; non-height-field loops and too-fine detail are refused. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01J53sR6gEaVt56soFDpTck8 * Test the Patch dialog filling a non-planar loop of picked edges The Patch dialog's only test used a sketch profile. Picking the four curved rim edges of a box cut by a sideways cylinder now makes one surface.patch command whose result is a single smooth face over the rim. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01J53sR6gEaVt56soFDpTck8 --------- Co-authored-by: Claude Opus 5.5 --- crates/engine/src/cmd/surface.rs | 23 +- crates/geom/src/bspline.rs | 178 +++--- crates/geom/src/fill.rs | 742 ++++++++++++++++++++++++++ crates/geom/src/lib.rs | 5 +- crates/kernel/src/fill.rs | 527 ++++++++++++++++++ crates/kernel/src/lib.rs | 3 + crates/kernel/src/plane_fit.rs | 133 +++++ crates/kernel/src/spline_conv.rs | 36 ++ crates/kernel/src/step_in/brep.rs | 28 +- crates/kernel/src/surfaces.rs | 74 +-- crates/kernel/src/tspline.rs | 20 +- crates/ui-egui/src/preselect_tests.rs | 25 + 12 files changed, 1648 insertions(+), 146 deletions(-) create mode 100644 crates/geom/src/fill.rs create mode 100644 crates/kernel/src/fill.rs create mode 100644 crates/kernel/src/plane_fit.rs create mode 100644 crates/kernel/src/spline_conv.rs diff --git a/crates/engine/src/cmd/surface.rs b/crates/engine/src/cmd/surface.rs index 0f430c90..0ee7b23e 100644 --- a/crates/engine/src/cmd/surface.rs +++ b/crates/engine/src/cmd/surface.rs @@ -13,7 +13,7 @@ pub static COMMANDS: &[CommandSpec] = &[ CommandSpec::new("surface.patch", "Patch", patch) .at("SURFACE", "CREATE") .icon("patch") - .params("sketch?: id|name and profiles? (planar patches of sketch regions) | edges: [[x,y,z]] (points on a closed loop of a body's edges) and body?; name?, body_name?"), + .params("sketch?: id|name and profiles? (planar patches of sketch regions) | edges: [[x,y,z]] (points on a closed loop of a body's edges: one face on them, planar when the loop is, else a smooth B-spline surface) and body?; name?, body_name?"), CommandSpec::new("surface.stitch", "Stitch", stitch) .at("SURFACE", "MODIFY") .icon("stitch") @@ -129,7 +129,8 @@ mod tests { solidcraft_kernel::measure(&st.body(name).unwrap().body).unwrap().volume } - /// A rectangle patched, trimmed by a plane and thickened; a box's faces stitched back. + /// A rectangle patched, trimmed by a plane and thickened; a box's open top patched; a rim + /// that is not planar patched. #[test] fn patch_trim_thicken_and_stitch() { let mut s = Session::default(); @@ -148,6 +149,24 @@ mod tests { run(&mut s, "solid.box", json!({"length": 10, "width": 10, "height": 10, "corner": [20, 0, 0], "body_name": "Cup"})); run(&mut s, "surface.patch", json!({"edges": [[25, 0, 10], [30, 5, 10], [25, 10, 10], [20, 5, 10]], "body": "Cup", "body_name": "Lid"})); assert!(s.world_state().body("Lid").unwrap().body.is_surface()); + // A rim that is not planar: a box with its top cut away by a sideways cylinder. + run(&mut s, "solid.box", json!({"length": 28, "width": 28, "height": 20, "corner": [-14, -14, 0], "body_name": "Saddle"})); + run( + &mut s, + "solid.cylinder", + json!({"radius": 40, "height": 120, "base": [0, -60, 50], "axis": "Y", "operation": "cut", "targets": ["Saddle"]}), + ); + let z = |x: f64| 50.0 - (1600.0 - x * x).sqrt(); + let rim = json!([[14, 0, z(14.0)], [0, 14, z(0.0)], [-14, 0, z(-14.0)], [0, -14, z(0.0)]]); + run(&mut s, "surface.patch", json!({"edges": rim, "body": "Saddle", "body_name": "Cover"})); + let st = s.world_state(); + let cover = &st.body("Cover").unwrap().body; + assert!(cover.is_surface() && cover.face_count() == 1); + // The cover is a smooth surface on the rim, close to the cylinder that cut it (its area + // within 1% of the cylinder's inside the rim). + let area = solidcraft_kernel::measure(cover).unwrap().area; + let want = 28.0 * 40.0 * 2.0 * (14.0f64 / 40.0).asin(); + assert!((area - want).abs() < want * 1e-2, "{area} vs {want}"); assert!(s.execute("surface.thicken", &json!({"bodies": ["Nope"], "thickness": 1})).is_err()); assert!(s.execute("surface.trim", &json!({"body": "Sheet"})).is_err()); } diff --git a/crates/geom/src/bspline.rs b/crates/geom/src/bspline.rs index a8e0f968..252a695e 100644 --- a/crates/geom/src/bspline.rs +++ b/crates/geom/src/bspline.rs @@ -154,6 +154,18 @@ pub fn check_knots(knots: &[f64]) -> SplineResult<()> { /// A knot vector shorter than two knots gives zeros. pub fn local_basis(knots: &[f64], u: f64, side: Side, n: usize) -> Vec { let mut out = vec![0.0; n + 1]; + let ders = local_basis_array(knots, u, side, n); + for (o, d) in out.iter_mut().zip(ders) { + *o = d; + } + out +} + +/// [`local_basis`] into a fixed array (no allocation): element `k` is the `k`-th derivative, +/// for `k ≤ min(n, MAX_SPLINE_DEGREE)`, the rest zero. +pub(crate) fn local_basis_array(knots: &[f64], u: f64, side: Side, n: usize) -> [f64; MAX_SPLINE_DEGREE + 1] { + const W: usize = MAX_SPLINE_DEGREE + 1; + let mut out = [0.0; W]; let Some(p) = knots.len().checked_sub(2) else { return out }; if p > MAX_SPLINE_DEGREE || !u.is_finite() { return out; @@ -163,7 +175,6 @@ pub fn local_basis(knots: &[f64], u: f64, side: Side, n: usize) -> Vec { // table[k][j]: k-th derivative of N_{j,q} for the current degree q (j = 0..=p−q), updated // in place: k descending (level k reads level k−1 of the previous degree), j ascending // (N_{j,q} reads N_{j+1,q−1}, not yet overwritten). Derivatives above p are zero. - const W: usize = MAX_SPLINE_DEGREE + 1; let levels = n.min(p); let mut table = [[0.0f64; W]; W]; for (j, slot) in table[0].iter_mut().enumerate().take(p + 1) { @@ -286,15 +297,16 @@ pub fn find_span(knots: &[f64], p: usize, n: usize, u: f64) -> Option { if u < lo || u > hi { return None; } - // The last non-empty span in [p, n−1] that starts at or before u. - let mut k = p; - for i in p..n { - if let (Some(&a), Some(&b)) = (knots.get(i), knots.get(i + 1)) - && a <= u - && a < b - { - k = i; - } + // The last non-empty span in [p, n−1] that starts at or before u: by bisection (the knots + // are sorted), and at the domain end the last non-empty span. + let mut k = if u >= hi { + n - 1 + } else { + let below = knots.get(p..n)?.partition_point(|&x| x <= u); + p + below.checked_sub(1)? + }; + while k > p && !(knots.get(k)? < knots.get(k + 1)?) { + k -= 1; } Some(k) } @@ -331,16 +343,62 @@ pub fn insert_knot(p: usize, knots: &[f64], ctrl: &[HPoint], u: f64) -> SplineRe Ok((nk, q)) } +/// The binomial coefficient `(n k)` (zero for `k > n`). +pub(crate) fn binomial(n: usize, k: usize) -> f64 { + if k > n { + return 0.0; + } + (0..k).fold(1.0, |c, m| c * (n - m) as f64 / (m + 1) as f64) +} + /// Bernstein polynomial `Bᵢ,ₚ(x)`. pub fn bernstein(p: usize, i: usize, x: f64) -> f64 { if i > p { return 0.0; } - let mut c = 1.0; - for m in 0..i { - c = c * (p - m) as f64 / (m + 1) as f64; + binomial(p, i) * x.powi(i as i32) * (1.0 - x).powi((p - i) as i32) +} + +/// One row of a B-spline (degree `p`, knots, homogeneous control points) clamped on its domain +/// `[knots[p], knots[n]]`: each domain end inserted until it repeats `p` times (§5.2), then the +/// basis functions that vanish on the domain dropped (a basis function's last span does not +/// depend on its first knot, so that knot moves to the domain end). The curve on the domain is +/// unchanged. +pub(crate) fn clamp_row(p: usize, knots: &[f64], row: &[HPoint]) -> SplineResult<(Vec, Vec)> { + let n = row.len(); + let (Some(&a), Some(&b)) = (knots.get(p), knots.get(n)) else { return invalid("too few knots") }; + let (mut knots, mut row) = (knots.to_vec(), row.to_vec()); + for (at_end, x) in [(false, a), (true, b)] { + let mult = |k: &[f64]| k.iter().filter(|&&y| y == x).count(); + let mut guard = 0; + while mult(&knots) < p { + guard += 1; + if guard > MAX_SPLINE_DEGREE + 1 { + return invalid("clamping did not converge"); + } + (knots, row) = insert_knot(p, &knots, &row, x)?; + } + let n = row.len(); + let (first, last) = if at_end { + // Functions i with i ≥ (first index of x) start at or past the end. + (0, knots.iter().position(|&y| y == x).unwrap_or(n).min(n)) + } else { + // Functions i with i + p + 1 ≤ (last index of x) end at or before the start. + (knots.iter().rposition(|&y| y == x).unwrap_or(p).saturating_sub(p), n) + }; + if last <= first + p { + return invalid("too few control points remain after clamping"); + } + let mut nk: Vec = knots.get(first..last + p + 1).map(<[f64]>::to_vec).unwrap_or_default(); + if at_end { + nk.last_mut().into_iter().for_each(|y| *y = x); + } else { + nk.first_mut().into_iter().for_each(|y| *y = x); + } + row = row.get(first..last).map(<[HPoint]>::to_vec).unwrap_or_default(); + knots = nk; } - c * x.powi(i as i32) * (1.0 - x).powi((p - i) as i32) + Ok((knots, row)) } /// A Bézier curve's point at `x ∈ [0, 1]` (de Casteljau, NURBS Book §1.4). @@ -543,66 +601,37 @@ impl NurbsSurface { } /// The same surface on the same domain with clamped knot vectors (`p+1` equal knots at - /// each end, so the boundary curves are given by the boundary rows of control points): - /// insert each domain end until it has multiplicity `p`, then drop the basis functions - /// that vanish on the domain (NURBS Book §5.2 knot insertion; a basis function's last span - /// does not depend on its first knot, so that knot can be moved to the domain end). + /// each end, so the boundary curves are given by the boundary rows of control points): each + /// row of control points clamped in turn ([`clamp_row`]). pub fn clamped(&self) -> SplineResult { let mut s = self.clone(); - for dir in [Dir::S, Dir::T] { - let d = dir.pick(0, 1); + for d in 0..2 { let p = s.degree[d]; - for end in 0..2 { - let a = s.domain()[d][end]; - while s.knots[d].iter().filter(|&&x| x == a).count() < p { - s = s.insert_knot(dir, a)?; - } - s = s.trim_end(d, end, a)?; + let [nu, nv] = s.count; + let (lines, len) = if d == 0 { (nv, nu) } else { (nu, nv) }; + let mut knots = Vec::new(); + let mut rows = Vec::with_capacity(lines); + for line in 0..lines { + let row: Vec = (0..len).map(|k| if d == 0 { s.point(k, line) } else { s.point(line, k) }).collect(); + let (k, r) = clamp_row(p, &s.knots[d], &row)?; + knots = k; + rows.push(r); } + let new_len = rows.first().map_or(0, Vec::len); + let mut ctrl = Vec::with_capacity(lines * new_len); + if d == 0 { + rows.iter().for_each(|r| ctrl.extend_from_slice(r)); + } else { + for j in 0..new_len { + ctrl.extend(rows.iter().map(|r| r.get(j).copied().unwrap_or(HPoint::ZERO))); + } + } + let mut all = s.knots.clone(); + all[d] = knots; + s = NurbsSurface::from_homogeneous(s.degree, all, ctrl)?; } Ok(s) } - - /// Drop the basis functions that vanish on the domain past the knot `a` (multiplicity ≥ p) - /// at the start (`end == 0`) or the end of direction `d`. - fn trim_end(&self, d: usize, end: usize, a: f64) -> SplineResult { - let p = self.degree[d]; - let kn = &self.knots[d]; - let n = self.count[d]; - let (first, last) = match end { - 0 => { - let last_a = kn.iter().rposition(|&x| x == a).ok_or_else(|| SplineError::Invalid("clamp".into()))?; - // Functions i with i + p + 1 ≤ last_a vanish on [a, …). - (last_a.saturating_sub(p), n) - } - _ => { - let first_a = kn.iter().position(|&x| x == a).ok_or_else(|| SplineError::Invalid("clamp".into()))?; - // Functions i with i ≥ first_a vanish on (…, a]. - (0, first_a.min(n)) - } - }; - if last <= first + p { - return invalid("clamp: too few control points remain"); - } - let mut knots = self.knots.clone(); - let mut nk: Vec = kn.get(first..last + p + 1).map(<[f64]>::to_vec).unwrap_or_default(); - match end { - 0 => nk.first_mut().into_iter().for_each(|x| *x = a), - _ => nk.last_mut().into_iter().for_each(|x| *x = a), - } - knots[d] = nk; - let [nu, nv] = self.count; - let mut ctrl = Vec::new(); - for j in 0..nv { - for i in 0..nu { - let k = if d == 0 { i } else { j }; - if k >= first && k < last { - ctrl.push(self.point(i, j)); - } - } - } - NurbsSurface::from_homogeneous(self.degree, knots, ctrl) - } } #[cfg(test)] @@ -666,6 +695,21 @@ pub(crate) mod tests { assert!(worst < 1e-5, "derivative error {worst}"); } + #[test] + fn span_search_finds_the_last_non_empty_span_at_or_before_u() { + let p = 3; + let knots = [0.0, 0.0, 0.0, 0.0, 0.3, 0.3, 0.5, 0.9, 0.9, 0.9, 1.0, 1.0, 1.0, 1.0]; + let n = knots.len() - p - 1; + for k in 0..=1000 { + let u = k as f64 / 1000.0; + // The plain search: the last span [knots[i], knots[i+1]) with knots[i] ≤ u. + let want = (p..n).filter(|&i| knots[i] <= u && knots[i] < knots[i + 1]).max().unwrap(); + assert_eq!(find_span(&knots, p, n, u), Some(want), "u = {u}"); + } + assert_eq!(find_span(&knots, p, n, 1.5), None); + assert_eq!(find_span(&knots, p, n, f64::NAN), None); + } + #[test] fn single_function_knot_insertion_is_exact() { let mut rng = Lcg(7); diff --git a/crates/geom/src/fill.rs b/crates/geom/src/fill.rs new file mode 100644 index 00000000..bfde4ac1 --- /dev/null +++ b/crates/geom/src/fill.rs @@ -0,0 +1,742 @@ +//! Filling a closed loop of curves with one surface: a non-rational bicubic B-spline height +//! field over the loop's best-fit plane, fitted to the loop by least squares with a smoothness +//! term. A loop that cannot be seen as the boundary of a height field over that plane (it +//! crosses itself, folds back or runs along the plane's normal when seen along that normal) is +//! refused, never filled with wrong geometry. +//! +//! References: Piegl & Tiller, *The NURBS Book* (2nd ed.) — §9.4.1 and §9.4.3 (least squares +//! approximation: minimise the sum of squared residuals, eq. 9.62, through the normal +//! equations, eq. 9.65), §3.4 (tensor product B-spline surfaces), §3.3 (derivative of a +//! B-spline curve, eq. 3.4: control values at the knot averages reproduce a linear function), +//! A2.1 (span search). + +use crate::Vec3; +use crate::bspline::{HPoint, NurbsSurface, Side, SplineError, SplineResult, find_span, invalid, local_basis_array}; + +/// Most loop points a fill takes. +pub const MAX_LOOP_POINTS: usize = 8192; +/// Most spans of the height field in each direction. +pub const MAX_FILL_SPANS: usize = 32; +/// Fewest spans of the height field in each direction (the first fit tried). +const MIN_FILL_SPANS: usize = 4; +/// The smoothness energy's terms: `(order in u, order in v, weight)` of `∫∫ (∂ h)²`. +const ENERGY: [(usize, usize, f64); 3] = [(2, 0, 1.0), (1, 1, 2.0), (0, 2, 1.0)]; +/// Four-point Gauss–Legendre nodes and weights on `[−1, 1]` (exact to degree 7). +const GAUSS4: [(f64, f64); 4] = [ + (-0.861_136_311_594_052_6, 0.347_854_845_137_453_8), + (-0.339_981_043_584_856_3, 0.652_145_154_862_546_2), + (0.339_981_043_584_856_3, 0.652_145_154_862_546_2), + (0.861_136_311_594_052_6, 0.347_854_845_137_453_8), +]; +/// The smoothness term's weight per fitted point (see [`fill_loop`]). +const SMOOTHNESS: f64 = 1e-7; +/// A loop step shorter than this fraction of its length in space, seen along the normal, runs +/// along the normal (steeper than about 87° from the plane). +const MIN_STEP_SLOPE: f64 = 0.05; +/// Consecutive steps closer to reversed than this cosine fold the loop back on itself. +const FOLD_COS: f64 = -0.999; + +/// A loop's fill: the surface (`S(u, v) = o + u·x + v·y + h(u, v)·n` over the plane with normal +/// `n`, so `∂S/∂u × ∂S/∂v` points along `n`), the plane's frame, whether the loop runs counter- +/// clockwise about it, the largest distance from a sampled loop point to the surface along `n`, +/// and how many points of the loop were evaluated. +#[derive(Clone, Debug, PartialEq)] +pub struct LoopFill { + pub surface: NurbsSurface, + /// The plane's origin and its axes `x`, `y` (the surface's parameter directions). + pub origin: Vec3, + pub axes: [Vec3; 2], + pub normal: Vec3, + pub counter_clockwise: bool, + pub deviation: f64, + pub evaluations: usize, +} + +/// Fill the closed loop of `sides` curves (side `k` given by `point(k, t)`, `t` from 0 at its +/// start to 1 at its end, which is where side `k+1` starts) with one surface passing within +/// `tol` of the loop. +/// +/// - Sampling: each side at `t = i/s` (fitted) and halfway between (checked), `s` from 16 +/// doubling while the fit needs denser points, at most [`MAX_LOOP_POINTS`] points in all. The +/// sample parameters are the same run either way, so the fill does not depend on where the +/// loop starts or which way it runs. +/// - The plane: through the fitted points' centroid, normal by Newell's method, turned so that +/// its largest component is positive; its first axis points to the fitted point farthest from +/// the centroid in the plane (the highest of points as far). So the plane and its axes, and with them the grid and the fill, +/// move with the loop (a reflected axis gives the same grid). +/// - Refused ("not supported yet"): a loop that, seen along the normal, has a step much shorter +/// than in space (it runs along the normal), turns straight back, or crosses or touches +/// itself; such a loop does not bound a height field over the plane. +/// - The height field: a bicubic B-spline with uniform knots over the projected loop's box +/// (widened by 1%); its control points sit at the knot averages in the plane, so its surface +/// parameters are plane coordinates, and their heights are the unknowns. The heights minimise +/// `Σ (h(uₖ, vₖ) − hₖ)² + λ·A · ∫∫ h_uu² + 2·h_uv² + h_vv²` over the fitted points — least +/// squares (eq. 9.62, normal equations eq. 9.65) plus the thin-plate bending energy over the +/// box, `A` the area of one grid cell (so each cell's energy counts like a squared second +/// difference of heights, and the balance is the same at any size). The energy decides the +/// inside, where there are no points; it is zero only for planes, so a plane through the +/// loop is reproduced exactly, and any other rim gets the least-bending surface on it, which +/// is free to turn at the rim (it is not the surface that made the rim). It does not depend +/// on the grid's directions. `λ = 10⁻⁷ ·` (fitted points) `/` (control points): just enough +/// to fix the inside, small enough not to pull the surface off the loop. +/// - Grids of 4, 8, 16, 32 spans across the box's longer side are tried, with as many across +/// the shorter side as keep the cells about square; a grid whose cells are under twice the +/// fitted points' largest step in a direction (cells between points would be left to the +/// energy alone) waits for denser sampling. The first grid within `tol` of every sampled +/// point (measured along the normal, which bounds the distance) is the fill; otherwise the +/// sampling doubles, and when that runs out, "not supported yet" with the closest deviation +/// reached. +pub fn fill_loop(sides: usize, point: &mut dyn FnMut(usize, f64) -> Option, tol: f64) -> SplineResult { + if sides == 0 || sides.saturating_mul(16) > MAX_LOOP_POINTS / 2 { + return Err(SplineError::TooLarge(format!("a loop of {sides} sides (at most {})", MAX_LOOP_POINTS / 32))); + } + if !(tol.is_finite() && tol > 0.0) { + return invalid("the fill tolerance must be positive"); + } + let mut evaluations = 0usize; + let mut eval = |k: usize, t: f64| -> SplineResult { + evaluations += 1; + point(k, t).filter(|p| p.is_finite()).ok_or_else(|| SplineError::Invalid(format!("side {k} has no finite point at {t}"))) + }; + let mut per = 16usize; + // Fitted points: side by side, per side at i/per (i < per). + let mut fit: Vec = Vec::with_capacity(sides * per); + for k in 0..sides { + for i in 0..per { + fit.push(eval(k, i as f64 / per as f64)?); + } + } + let mut deviation = f64::INFINITY; + loop { + let mut check = Vec::with_capacity(fit.len()); + for k in 0..sides { + for i in 0..per { + check.push(eval(k, (i as f64 + 0.5) / per as f64)?); + } + } + // The loop in order: fitted and checked points alternating. + let all: Vec = fit.iter().zip(&check).flat_map(|(a, b)| [*a, *b]).collect(); + match fit_once(&fit, &all, tol)? { + Fit::Done(mut fill) => { + fill.evaluations = evaluations; + return Ok(*fill); + } + Fit::Denser(dev) => deviation = deviation.min(dev), + Fit::Exhausted(dev) => { + deviation = deviation.min(dev); + break; + } + } + if all.len() * 2 > MAX_LOOP_POINTS { + break; + } + fit = all; + per *= 2; + } + Err(SplineError::Unsupported(format!("filling a loop this wavy with one surface (it stays {deviation:e} off the loop, tolerance {tol:e})"))) +} + +/// What one round of [`fill_loop`] gave. +enum Fit { + Done(Box), + /// Within tolerance only with cells finer than the points allow: sample denser. + Denser(f64), + /// The finest grid tried is not within tolerance. + Exhausted(f64), +} + +/// One round of [`fill_loop`]: fit the points `fit`, check every point of `all` (the loop in +/// order, `fit` among them). +fn fit_once(fit: &[Vec3], all: &[Vec3], tol: f64) -> SplineResult { + let m = fit.len(); + let o = fit.iter().fold(Vec3::ZERO, |a, p| a + *p) / m as f64; + let at = |k: usize| fit.get(k % m).copied().unwrap_or(Vec3::ZERO); + let area = (0..m).fold(Vec3::ZERO, |a, k| a + (at(k) - o).cross(at(k + 1) - o)); + let Some(normal) = area.normalized() else { return invalid("the loop encloses no area") }; + let largest = [normal.x, normal.y, normal.z].into_iter().fold(0.0f64, |a, c| if c.abs() > a.abs() { c } else { a }); + let counter_clockwise = largest > 0.0; + let n = if counter_clockwise { normal } else { normal * -1.0 }; + // The first axis towards the fitted point farthest from the centroid in the plane: the axes, + // the grid and so the fill move with the loop. + let flat = |p: &Vec3| { + let d = *p - o; + d - n * d.dot(n) + }; + // Points as far as each other (within rounding) are told apart by their height, and points + // alike in both by their coordinates: the same point is chosen however the loop is given, + // and as the loop moves, but for ties its symmetry makes. + let key = |p: &Vec3| (flat(p).len(), (*p - o).dot(n)); + let reach = fit.iter().map(|p| key(p).0).fold(0.0, f64::max); + let near = 1e-9 * reach; + let far = fit.iter().fold(None::<&Vec3>, |best, p| match best { + None => Some(p), + Some(b) => { + let ((db, hb), (dp, hp)) = (key(b), key(p)); + let better = if (dp - db).abs() > near { + dp > db + } else if (hp - hb).abs() > near { + hp > hb + } else { + (p.x, p.y, p.z) > (b.x, b.y, b.z) + }; + Some(if better { p } else { b }) + } + }); + let far = far.map(flat).unwrap_or(Vec3::ZERO); + let x = far.normalized().ok_or_else(|| SplineError::Invalid("the loop has no extent in its plane".into()))?; + let y = n.cross(x); + let project = |pts: &[Vec3]| -> (Vec<[f64; 2]>, Vec) { + (pts.iter().map(|p| [(*p - o).dot(x), (*p - o).dot(y)]).collect(), pts.iter().map(|p| (*p - o).dot(n)).collect()) + }; + let (uv, h) = project(all); + check_projectable(all, &uv)?; + let (fuv, fh) = project(fit); + let (mut lo, mut hi) = ([f64::INFINITY; 2], [f64::NEG_INFINITY; 2]); + for q in &uv { + for d in 0..2 { + lo[d] = lo[d].min(q[d]); + hi[d] = hi[d].max(q[d]); + } + } + for d in 0..2 { + let pad = 0.01 * (hi[d] - lo[d]); + lo[d] -= pad; + hi[d] += pad; + } + if !(lo[0] < hi[0] && lo[1] < hi[1] && (hi[0] - lo[0]).is_finite() && (hi[1] - lo[1]).is_finite()) { + return invalid("the loop's extent is not finite"); + } + // Spans per direction: `level` across the longer side of the box, as many across the + // shorter as keep the cells about square (at least one). A grid finer than the fitted + // points can hold in a direction (cells under twice their largest step along it) is tried + // only after denser sampling; the coarsest is always tried. + let extent = [hi[0] - lo[0], hi[1] - lo[1]]; + let mut step = [0.0f64; 2]; + for k in 0..m { + let (p, q) = (fuv.get(k).copied().unwrap_or([0.0; 2]), fuv.get((k + 1) % m).copied().unwrap_or([0.0; 2])); + for d in 0..2 { + step[d] = step[d].max((q[d] - p[d]).abs()); + } + } + let longer = extent[0].max(extent[1]); + let mut deviation = f64::INFINITY; + let mut level = MIN_FILL_SPANS; + while level <= MAX_FILL_SPANS { + let spans = [0, 1].map(|d| ((extent[d] / longer * level as f64).ceil() as usize).clamp(1, level)); + let holds = (0..2).all(|d| extent[d] / spans[d] as f64 >= 2.0 * step[d]); + if !holds && level > MIN_FILL_SPANS { + return Ok(Fit::Denser(deviation)); + } + let grid = Grid::new(spans, lo, hi); + let heights = grid.fit(&fuv, &fh)?; + let dev = uv.iter().zip(&h).map(|(q, hk)| (grid.height(&heights, *q) - hk).abs()).fold(0.0, f64::max); + if !dev.is_finite() { + return invalid("the fit is not finite"); + } + deviation = deviation.min(dev); + if dev <= tol && holds { + let surface = grid.surface(&heights, o, x, y, n)?; + let fill = LoopFill { surface, origin: o, axes: [x, y], normal: n, counter_clockwise, deviation: dev, evaluations: 0 }; + return Ok(Fit::Done(Box::new(fill))); + } + if !holds { + return Ok(Fit::Denser(deviation)); + } + level *= 2; + } + Ok(Fit::Exhausted(deviation)) +} + +/// Refuse a loop that does not bound a height field over its plane (see [`fill_loop`]). +fn check_projectable(points: &[Vec3], uv: &[[f64; 2]]) -> SplineResult<()> { + let m = uv.len(); + let q = |k: usize| uv.get(k % m).copied().unwrap_or([0.0; 2]); + let p = |k: usize| points.get(k % m).copied().unwrap_or(Vec3::ZERO); + let step = |k: usize| { + let (a, b) = (q(k), q(k + 1)); + [b[0] - a[0], b[1] - a[1]] + }; + for k in 0..m { + let d = step(k); + let flat = d[0].hypot(d[1]); + let space = p(k).dist(p(k + 1)); + if space > 0.0 && flat < MIN_STEP_SLOPE * space { + return Err(SplineError::Unsupported("filling a loop that runs along its plane's normal".into())); + } + let e = step(k + 1); + let fe = e[0].hypot(e[1]); + if flat > 0.0 && fe > 0.0 && (d[0] * e[0] + d[1] * e[1]) < FOLD_COS * flat * fe { + return Err(SplineError::Unsupported("filling a loop that folds back on itself".into())); + } + } + // Non-adjacent segments must not meet: a sweep over the segments' u extents. + let mut segs: Vec<(f64, f64, usize)> = (0..m) + .map(|k| { + let (a, b) = (q(k), q(k + 1)); + (a[0].min(b[0]), a[0].max(b[0]), k) + }) + .collect(); + segs.sort_by(|a, b| a.0.total_cmp(&b.0)); + for (i, &(_, hi, k)) in segs.iter().enumerate() { + for &(lo2, _, l) in segs.get(i + 1..).into_iter().flatten() { + if lo2 > hi { + break; + } + let adjacent = (k + 1) % m == l || (l + 1) % m == k; + if !adjacent && segments_meet(q(k), q(k + 1), q(l), q(l + 1)) { + return Err(SplineError::Unsupported("filling a loop that crosses itself seen along its plane's normal".into())); + } + } + } + Ok(()) +} + +/// Whether the closed segments `ab` and `cd` share a point. +fn segments_meet(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> bool { + let orient = |p: [f64; 2], q: [f64; 2], r: [f64; 2]| (q[0] - p[0]) * (r[1] - p[1]) - (q[1] - p[1]) * (r[0] - p[0]); + let within = + |p: [f64; 2], q: [f64; 2], r: [f64; 2]| r[0] >= p[0].min(q[0]) && r[0] <= p[0].max(q[0]) && r[1] >= p[1].min(q[1]) && r[1] <= p[1].max(q[1]); + let (d1, d2, d3, d4) = (orient(c, d, a), orient(c, d, b), orient(a, b, c), orient(a, b, d)); + if ((d1 > 0.0 && d2 < 0.0) || (d1 < 0.0 && d2 > 0.0)) && ((d3 > 0.0 && d4 < 0.0) || (d3 < 0.0 && d4 > 0.0)) { + return true; + } + (d1 == 0.0 && within(c, d, a)) || (d2 == 0.0 && within(c, d, b)) || (d3 == 0.0 && within(a, b, c)) || (d4 == 0.0 && within(a, b, d)) +} + +/// A bicubic height field's grid: uniform clamped knots over `[lo, hi]` in each direction, +/// `count` control points along each. +struct Grid { + knots: [Vec; 2], + count: [usize; 2], +} + +impl LoopFill { + /// The surface's point over (or under) `p`, along the normal. + pub fn point_over(&self, p: Vec3) -> SplineResult { + let d = p - self.origin; + let [x, y] = self.axes; + self.surface.eval(d.dot(x), d.dot(y)) + } +} + +impl Grid { + fn new(spans: [usize; 2], lo: [f64; 2], hi: [f64; 2]) -> Grid { + let knots = [0, 1].map(|d| { + let (a, b, s) = (lo[d], hi[d], spans[d]); + let mut k = vec![a; 4]; + k.extend((1..s).map(|i| a + (b - a) * i as f64 / s as f64)); + k.extend([b; 4]); + k + }); + Grid { knots, count: spans.map(|s| s + 3) } + } + + /// The length of one span in direction `d`. + fn cell(&self, d: usize) -> f64 { + let (Some(k), Some(&c)) = (self.knots.get(d), self.count.get(d)) else { return 0.0 }; + (k.last().copied().unwrap_or(0.0) - k.first().copied().unwrap_or(0.0)) / (c - 3) as f64 + } + + /// The knot average `(kᵢ₊₁ + kᵢ₊₂ + kᵢ₊₃)/3` of control point `i` in direction `d`. + fn average(&self, d: usize, i: usize) -> f64 { + self.knots.get(d).map_or(0.0, |k| (1..=3).map(|t| k.get(i + t).copied().unwrap_or(0.0)).sum::() / 3.0) + } + + /// The four non-zero basis functions' `r`-th derivatives and the first index, at `x` in + /// direction `d`. + fn basis_derivative(&self, d: usize, x: f64, r: usize) -> Option<(usize, [f64; 4])> { + let kn = self.knots.get(d)?; + let x = x.clamp(*kn.first()?, *kn.last()?); + let span = find_span(kn, 3, *self.count.get(d)?, x)?; + let side = if x >= *kn.last()? { Side::Left } else { Side::Right }; + let first = span.checked_sub(3)?; + let mut out = [0.0; 4]; + for (a, o) in out.iter_mut().enumerate() { + *o = kn.get(first + a..first + a + 5).map_or(0.0, |l| local_basis_array(l, x, side, r).get(r).copied().unwrap_or(0.0)); + } + Some((first, out)) + } + + /// The height at `q` of the net `heights` (i fastest). + fn height(&self, heights: &[f64], q: [f64; 2]) -> f64 { + let (Some((i0, bu)), Some((j0, bv))) = (self.basis_derivative(0, q[0], 0), self.basis_derivative(1, q[1], 0)) else { return f64::NAN }; + let mut s = 0.0; + for (b, wv) in bv.iter().enumerate() { + for (a, wu) in bu.iter().enumerate() { + s += wu * wv * heights.get((j0 + b) * self.count[0] + i0 + a).copied().unwrap_or(f64::NAN); + } + } + s + } + + /// The control heights fitting the points (see [`fill_loop`]): the normal equations, banded + /// (unknowns i fastest, so coupled unknowns are at most `3·count[0] + 3` apart), solved by + /// Cholesky factorisation. + fn fit(&self, uv: &[[f64; 2]], h: &[f64]) -> SplineResult> { + let [cu, cv] = self.count; + let n = cu * cv; + let mut band = Band::new(n, 3 * cu + 3); + let mut rhs = vec![0.0; n]; + let mut fitted = 0usize; + for (q, hk) in uv.iter().zip(h) { + let (Some((i0, bu)), Some((j0, bv))) = (self.basis_derivative(0, q[0], 0), self.basis_derivative(1, q[1], 0)) else { + return invalid("a loop point is off the grid"); + }; + let mut terms = [(0usize, 0.0f64); 16]; + for (t, slot) in terms.iter_mut().enumerate() { + let (a, b) = (t % 4, t / 4); + *slot = ((j0 + b) * cu + i0 + a, bu.get(a).copied().unwrap_or(0.0) * bv.get(b).copied().unwrap_or(0.0)); + } + for &(r, wr) in &terms { + if let Some(slot) = rhs.get_mut(r) { + *slot += wr * hk; + } + for &(s, ws) in &terms { + band.add(r, s, wr * ws); + } + } + fitted += 1; + } + // The thin-plate energy ∫∫ h_uu² + 2·h_uv² + h_vv² over the grid's box, times the area of + // one cell: per cell that is the squared second difference of the heights, so it weighs + // against the data term (squared heights) alike at any size. Each term is a product of + // one-direction Gram matrices of the basis functions' derivatives. + let lambda = SMOOTHNESS * fitted as f64 / n as f64 * self.cell(0) * self.cell(1); + let gram = |d: usize, r: usize| -> Vec { + let c = self.count.get(d).copied().unwrap_or(0); + let mut g = vec![0.0; c * c]; + let Some(k) = self.knots.get(d) else { return g }; + for w in k.windows(2) { + let &[a, b] = w else { continue }; + if !(b > a) { + continue; + } + for (gx, gw) in GAUSS4 { + let x = 0.5 * (a + b) + 0.5 * (b - a) * gx; + let Some((first, vals)) = self.basis_derivative(d, x, r) else { continue }; + for (i, vi) in vals.iter().enumerate() { + for (j, vj) in vals.iter().enumerate() { + if let Some(slot) = g.get_mut((first + i) * c + first + j) { + *slot += 0.5 * (b - a) * gw * vi * vj; + } + } + } + } + } + g + }; + let grams: [[Vec; 3]; 2] = [0, 1].map(|d| [0, 1, 2].map(|r| gram(d, r))); + let g = |d: usize, r: usize, i: usize, k: usize| { + let c = self.count.get(d).copied().unwrap_or(0); + grams.get(d).and_then(|m| m.get(r)).and_then(|m| m.get(i * c + k)).copied().unwrap_or(0.0) + }; + for j in 0..cv { + for i in 0..cu { + let r = j * cu + i; + for l in j.saturating_sub(3)..(j + 4).min(cv) { + for k in i.saturating_sub(3)..(i + 4).min(cu) { + let e: f64 = ENERGY.iter().map(|&(ru, rv, w)| w * g(0, ru, i, k) * g(1, rv, j, l)).sum(); + band.add(r, l * cu + k, lambda * e); + } + } + } + } + band.solve(rhs) + } + + /// The surface of the net `heights` over the plane `o`, `x`, `y`, `n`: control point `(i, j)` + /// at the knot averages `(ξᵢ, ηⱼ)` in the plane, `heights[j·count[0] + i]` along `n`. + fn surface(&self, heights: &[f64], o: Vec3, x: Vec3, y: Vec3, n: Vec3) -> SplineResult { + let [cu, cv] = self.count; + let mut ctrl = Vec::with_capacity(cu * cv); + for j in 0..cv { + for i in 0..cu { + let hij = heights.get(j * cu + i).copied().unwrap_or(f64::NAN); + ctrl.push(HPoint::weighted(o + x * self.average(0, i) + y * self.average(1, j) + n * hij, 1.0)); + } + } + NurbsSurface::from_homogeneous([3, 3], self.knots.clone(), ctrl) + } +} + +/// A symmetric banded matrix: the upper band, row `r` holding columns `r ..= r + bw`. +struct Band { + n: usize, + bw: usize, + a: Vec, +} + +impl Band { + fn new(n: usize, bw: usize) -> Band { + Band { n, bw, a: vec![0.0; n * (bw + 1)] } + } + fn index(&self, r: usize, s: usize) -> Option { + let (r, s) = if r <= s { (r, s) } else { (s, r) }; + (s - r <= self.bw && s < self.n).then_some(r * (self.bw + 1) + (s - r)) + } + fn get(&self, r: usize, s: usize) -> f64 { + self.index(r, s).and_then(|i| self.a.get(i)).copied().unwrap_or(0.0) + } + fn set(&mut self, r: usize, s: usize, v: f64) { + if let Some(x) = self.index(r, s).and_then(|i| self.a.get_mut(i)) { + *x = v; + } + } + /// Add `v` at `(r, s)` with `r ≤ s` (the lower half is the same by symmetry: callers adding + /// a full outer product reach each stored entry once through `r ≤ s`). + fn add(&mut self, r: usize, s: usize, v: f64) { + if r <= s + && let Some(x) = self.index(r, s).and_then(|i| self.a.get_mut(i)) + { + *x += v; + } + } + /// Solve `A·x = b` by banded Cholesky factorisation `A = Uᵀ·U`; an error when `A` is not + /// positive definite. + fn solve(mut self, mut b: Vec) -> SplineResult> { + let (n, bw) = (self.n, self.bw); + for r in 0..n { + let mut d = self.get(r, r); + for k in r.saturating_sub(bw)..r { + let u = self.get(k, r); + d -= u * u; + } + if !(d > 0.0) { + return invalid("the fit's equations are singular"); + } + let d = d.sqrt(); + self.set(r, r, d); + for s in r + 1..(r + bw + 1).min(n) { + let mut v = self.get(r, s); + for k in s.saturating_sub(bw)..r { + v -= self.get(k, r) * self.get(k, s); + } + self.set(r, s, v / d); + } + } + // Uᵀ·y = b, then U·x = y. + for r in 0..n { + let mut v = b.get(r).copied().unwrap_or(0.0); + for k in r.saturating_sub(bw)..r { + v -= self.get(k, r) * b.get(k).copied().unwrap_or(0.0); + } + if let Some(x) = b.get_mut(r) { + *x = v / self.get(r, r); + } + } + for r in (0..n).rev() { + let mut v = b.get(r).copied().unwrap_or(0.0); + for s in r + 1..(r + bw + 1).min(n) { + v -= self.get(r, s) * b.get(s).copied().unwrap_or(0.0); + } + if let Some(x) = b.get_mut(r) { + *x = v / self.get(r, r); + } + } + if b.iter().any(|x| !x.is_finite()) { + return invalid("the fit is not finite"); + } + Ok(b) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// The closed polygon through `corners` in space, side `k` from corner `k` to `k+1`. + fn polygon3(corners: &[Vec3]) -> impl FnMut(usize, f64) -> Option + '_ { + move |k, t| { + let n = corners.len(); + Some(corners.get(k)?.lerp(*corners.get((k + 1) % n)?, t)) + } + } + + /// The polygon `corners` in the plane, lifted by `z(x, y)` all along. + fn lifted<'a>(corners: &'a [[f64; 2]], z: &'a dyn Fn(f64, f64) -> f64) -> impl FnMut(usize, f64) -> Option + 'a { + move |k, t| { + let n = corners.len(); + let (a, b) = (corners.get(k)?, corners.get((k + 1) % n)?); + let (x, y) = (a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t); + Some(Vec3::new(x, y, z(x, y))) + } + } + + fn ring(n: usize, r: f64) -> Vec<[f64; 2]> { + (0..n).map(|k| k as f64 * std::f64::consts::TAU / n as f64).map(|a| [r * a.cos(), r * a.sin()]).collect() + } + + /// The farthest the surface is from `z` over the polygon's inside (sampled on a grid). + fn off_inside(f: &LoopFill, corners: &[[f64; 2]], z: &dyn Fn(f64, f64) -> f64) -> f64 { + let inside = |x: f64, y: f64| { + let n = corners.len(); + let mut w = false; + for k in 0..n { + let (a, b) = (corners[k], corners[(k + 1) % n]); + if (a[1] > y) != (b[1] > y) && x < a[0] + (y - a[1]) / (b[1] - a[1]) * (b[0] - a[0]) { + w = !w; + } + } + w + }; + let s = &f.surface; + let [du, dv] = s.domain(); + let mut worst = 0.0f64; + for i in 0..=40 { + for j in 0..=40 { + let (u, v) = ((du[0] + (du[1] - du[0]) * i as f64 / 40.0).min(du[1]), (dv[0] + (dv[1] - dv[0]) * j as f64 / 40.0).min(dv[1])); + let p = s.eval(u, v).unwrap(); + if inside(p.x, p.y) { + worst = worst.max((p.z - z(p.x, p.y)).abs()); + } + } + } + worst + } + + #[test] + fn planes_are_reproduced_and_rims_of_quadrics_are_filled_closely() { + // A tilted plane over a hexagon: the smoothness energy vanishes on it, so the fill is it. + let hex = ring(6, 10.0); + let plane = |x: f64, y: f64| 0.3 * x - 0.2 * y + 1.0; + let f = fill_loop(6, &mut lifted(&hex, &plane), 1e-6).unwrap(); + let inside = off_inside(&f, &hex, &plane); + eprintln!("plane: deviation {:e}, inside off by {inside:e}, {} evaluations", f.deviation, f.evaluations); + assert!(f.deviation <= 1e-6 && inside < 1e-9); + // A saddle and a dome: the least-bending fill through the rim is not the quadric (only + // planes are reproduced), but stays near it (the plate is free to turn at the rim). + let saddle = |x: f64, y: f64| 0.05 * x * y; + let dome = |x: f64, y: f64| 3.0 - 0.02 * (x * x + y * y); + for (name, z, bound) in [("saddle", &saddle as &dyn Fn(f64, f64) -> f64, 0.1), ("dome", &dome, 1.0)] { + let f = fill_loop(6, &mut lifted(&hex, z), 2e-5).unwrap(); + let inside = off_inside(&f, &hex, z); + eprintln!("{name}: deviation {:e}, inside off by {inside:e}", f.deviation); + assert!(f.deviation <= 2e-5 && inside < bound, "{name}: {inside}"); + } + } + + #[test] + fn concave_loops_fill_and_the_result_does_not_depend_on_where_the_loop_starts() { + // A dart, an arrow-head pentagon and a plus, on a cylinder. + let z = |x: f64, _: f64| 50.0 - (1600.0 - x * x).sqrt(); + let dart = vec![[0.0, 0.0], [20.0, -10.0], [8.0, 0.0], [20.0, 10.0]]; + let arrow = vec![[0.0, 0.0], [10.0, -8.0], [10.0, -3.0], [20.0, 0.0], [10.0, 3.0]]; + let plus = vec![ + [-3.0, -9.0], + [3.0, -9.0], + [3.0, -3.0], + [9.0, -3.0], + [9.0, 3.0], + [3.0, 3.0], + [3.0, 9.0], + [-3.0, 9.0], + [-3.0, 3.0], + [-9.0, 3.0], + [-9.0, -3.0], + [-3.0, -3.0], + ]; + let flat_corner = vec![[-10.0, -10.0], [0.0, -10.0], [10.0, -10.0], [10.0, 10.0], [-10.0, 10.0]]; + for (name, c) in [("dart", dart), ("arrow", arrow), ("plus", plus), ("flat corner", flat_corner)] { + let n = c.len(); + let f = fill_loop(n, &mut lifted(&c, &z), 1e-4).unwrap(); + let inside = off_inside(&f, &c, &z); + eprintln!("{name}: deviation {:e}, inside off the cylinder by {inside:e}, {} evaluations", f.deviation, f.evaluations); + // The least-bending fill, not the cylinder: within a few tenths of a millimetre of it. + assert!(f.deviation <= 1e-4 && inside < 0.2, "{name}: {inside}"); + // The same loop from another corner and the other way round. + let mut back = c.clone(); + back.rotate_left(2); + back.reverse(); + let g = fill_loop(n, &mut lifted(&back, &z), 1e-4).unwrap(); + assert_ne!(f.counter_clockwise, g.counter_clockwise); + for (x, y) in [(4.0, 0.0), (6.0, 1.0), (5.0, -1.0)] { + let p = Vec3::new(x, y, z(x, y)); + let gap = f.point_over(p).unwrap().dist(g.point_over(p).unwrap()); + assert!(gap < 1e-7, "{name}: {gap}"); + } + } + } + + #[test] + fn the_fill_is_the_same_at_any_size_and_place() { + // The arrow on the cylinder, scaled by 100 and by 1/100, turned and moved: the same fill, + // scaled and moved alike. + let z = |x: f64, _: f64| 50.0 - (1600.0 - x * x).sqrt(); + let arrow = [[0.0, 0.0], [10.0, -8.0], [10.0, -3.0], [20.0, 0.0], [10.0, 3.0]]; + let f = fill_loop(5, &mut lifted(&arrow, &z), 1e-4).unwrap(); + let axis = Vec3::new(1.0, -2.0, 0.5).normalized().unwrap(); + let shift = Vec3::new(30.0, -12.0, 7.0); + for scale in [100.0, 0.01] { + let place = |p: Vec3| { + // Scaled, turned by 0.9 about `axis` (Rodrigues), moved. + let p = p * scale; + let (c, s) = (0.9f64.cos(), 0.9f64.sin()); + p * c + axis.cross(p) * s + axis * (axis.dot(p) * (1.0 - c)) + shift * scale + }; + let mut moved = |k: usize, t: f64| lifted(&arrow, &z)(k, t).map(place); + let g = fill_loop(5, &mut moved, 1e-4 * scale).unwrap(); + let mut worst = 0.0f64; + for (x, y) in [(4.0, 0.0), (6.0, 1.0), (5.0, -1.0), (9.0, -2.0), (12.0, 0.5)] { + let p = Vec3::new(x, y, z(x, y)); + worst = worst.max(g.point_over(place(p)).unwrap().dist(place(f.point_over(p).unwrap())) / scale); + } + eprintln!("arrow scaled by {scale}, turned and moved: {} evaluations, surfaces apart by {worst:e} (relative)", g.evaluations); + assert!(worst < 1e-7, "scale {scale}: {worst}"); + } + } + + #[test] + fn a_long_narrow_slot_fills() { + // A 200 × 1 mm slot on a cylinder of radius 300 running along it. + let z = |x: f64, _: f64| 300.0 - (90_000.0 - (x - 100.0) * (x - 100.0)).sqrt(); + let slot = [[0.0, 0.0], [200.0, 0.0], [200.0, 1.0], [0.0, 1.0]]; + let f = fill_loop(4, &mut lifted(&slot, &z), 2e-4).unwrap(); + let inside = off_inside(&f, &slot, &z); + eprintln!("slot: deviation {:e}, inside off the cylinder by {inside:e}, {} evaluations", f.deviation, f.evaluations); + assert!(f.deviation <= 2e-4 && inside < 1e-2); + } + + #[test] + fn loops_that_do_not_bound_a_height_field_are_refused() { + let unsupported = |r: SplineResult| matches!(r, Err(SplineError::Unsupported(_))); + let flat = |_: f64, _: f64| 0.0; + // A figure eight (lobes of different sizes). + let eight = [[0.0, 0.0], [10.0, 10.0], [10.0, 0.0], [0.0, 6.0]]; + assert!(unsupported(fill_loop(4, &mut lifted(&eight, &flat), 1e-4))); + // A loop with a vertical spike in it (up along the normal and down again). + let spike = [ + Vec3::ZERO, + Vec3::new(10.0, 0.0, 0.0), + Vec3::new(10.0, 10.0, 0.0), + Vec3::new(5.0, 10.0, 0.0), + Vec3::new(5.0, 10.0, 3.0), + Vec3::new(5.0, 10.0, 0.0), + Vec3::new(0.0, 10.0, 0.0), + ]; + assert!(unsupported(fill_loop(7, &mut polygon3(&spike), 1e-4))); + // A loop that runs back over itself. + let back = [[0.0, 0.0], [10.0, 0.0], [4.0, 0.0], [10.0, 10.0], [0.0, 10.0]]; + assert!(unsupported(fill_loop(5, &mut lifted(&back, &flat), 1e-4))); + // A loop too wavy for one surface. + let wavy = |x: f64, y: f64| (x * 9.0).sin() * (y * 7.0).cos(); + assert!(unsupported(fill_loop(64, &mut lifted(&ring(64, 10.0), &wavy), 1e-6))); + } + + #[test] + fn hostile_loops_are_refused() { + let flat = |_: f64, _: f64| 0.0; + let five = ring(5, 10.0); + assert!(fill_loop(0, &mut lifted(&five, &flat), 1e-4).is_err()); + assert!(fill_loop(5, &mut |_, _| Some(Vec3::X), 1e-4).is_err(), "no area"); + assert!(fill_loop(5, &mut |_, _| None, 1e-4).is_err()); + assert!(fill_loop(5, &mut |k, t| Some(Vec3::new(k as f64, t, f64::NAN)), 1e-4).is_err()); + assert!(fill_loop(5, &mut lifted(&five, &flat), 0.0).is_err() && fill_loop(5, &mut lifted(&five, &flat), f64::NAN).is_err()); + assert!(fill_loop(10_000, &mut lifted(&five, &flat), 1e-4).is_err()); + let huge = |x: f64, y: f64| (x * y) * 1e300; + let _ = fill_loop(5, &mut lifted(&five, &huge), 1e-4); + let tiny: Vec<[f64; 2]> = five.iter().map(|p| [p[0] * 1e-300, p[1] * 1e-300]).collect(); + let _ = fill_loop(5, &mut lifted(&tiny, &flat), 1e-4); + let big: Vec<[f64; 2]> = five.iter().map(|p| [p[0] * 1e300, p[1] * 1e300]).collect(); + let _ = fill_loop(5, &mut lifted(&big, &flat), 1e-4); + } +} diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index ec9d519b..d09a07d0 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -1,7 +1,7 @@ //! Geometry foundation for SolidCraft: 2D/3D vectors, sketch planes, axis-aligned boxes, //! triangle meshes with their measures (volume, area, centroid), and freeform surface maths -//! (B-spline basis and NURBS surfaces in [`bspline`], T-splines in [`tspline`]). Units are -//! millimetres. +//! (B-spline basis and NURBS surfaces in [`bspline`], one-surface fills of closed loops in +//! [`fill`], T-splines in [`tspline`]). Units are millimetres. //! //! This crate has no knowledge of sketches, documents or the B-rep kernel; every other crate //! builds on it. @@ -10,6 +10,7 @@ pub mod bspline; pub mod env; +pub mod fill; mod mesh; mod plane; mod profile; diff --git a/crates/kernel/src/fill.rs b/crates/kernel/src/fill.rs new file mode 100644 index 00000000..829d84c8 --- /dev/null +++ b/crates/kernel/src/fill.rs @@ -0,0 +1,527 @@ +//! Filling a closed loop of a body's edges that is not planar with one face: a B-spline height +//! field over the loop's best-fit plane (`solidcraft_geom::fill::fill_loop`), bounded by the +//! body's own edges, whole. Loops that do not bound a height field over a plane are refused. + +use solidcraft_geom::Vec3; +use solidcraft_geom::bspline::SplineError; +use truck_modeling as mt; + +use crate::body::{Body, from_p3}; +use crate::plane_fit::EdgeSampler; +use crate::record::{Rebuild, rebuilt_body}; +use crate::{KernelError, Result}; + +/// A surface body of one face filling the closed loop `chain` of `body`'s edges (each edge +/// oriented along the loop, each ending where the next starts), within [`face_tolerance`] of +/// every edge. Its history is over `[body]`: the face is new, its edges the body's own. Call it +/// under [`crate::guard`]. +pub(crate) fn fill_loop(body: &Body, chain: &[mt::Edge]) -> Result { + let tol = face_tolerance(body); + let mut sides: Vec<(EdgeSampler, Vec3)> = Vec::with_capacity(chain.len()); + for e in chain { + let sampler = EdgeSampler::new(e).ok_or_else(|| KernelError::Invalid("an edge has no parameter range".into()))?; + sides.push((sampler, from_p3(e.front().point()))); + } + let mut point = |k: usize, t: f64| -> Option { + let (sampler, start) = sides.get_mut(k)?; + // The loop's corners are its vertices (an imported edge's curve may end a little off). + if t <= 0.0 { Some(*start) } else { sampler.at(t) } + }; + let fill = solidcraft_geom::fill::fill_loop(chain.len(), &mut point, tol).map_err(|e| match e { + SplineError::Unsupported(m) => KernelError::Invalid(format!("not supported yet: {m}")), + other => KernelError::Invalid(format!("cannot fill the loop: {other}")), + })?; + // The wire runs counter-clockwise about the surface's normal. + let wire: mt::Wire = + if fill.counter_clockwise { chain.to_vec().into() } else { chain.iter().rev().map(mt::Edge::inverse).collect::>().into() }; + let face = mt::Face::try_new(vec![wire], crate::spline_conv::surface(&fill.surface)?).map_err(|e| KernelError::Failed(format!("fill: {e}")))?; + let mut rec = Rebuild::default(); + rec.made(&face, Vec::new()); + Ok(rebuilt_body(crate::surfaces::surface_body(vec![face]), &[&body.solid], &rec)) +} + +/// How far a loop's edges may be from a plane for its patch to be that plane: `10⁻⁷` of the +/// body's size, at least the linear tolerance (`10⁻⁶` mm) — the body's own "same point" +/// tolerance, since a planar face is exact and other operations treat it as containing its +/// edges to that tolerance. A fitted face gets [`face_tolerance`], ten times looser. +pub(crate) fn planarity_tolerance(body: &Body) -> f64 { + (body.size() * 1e-7).max(solidcraft_geom::LINEAR_TOL) +} + +/// How far the fill's surface may be from the loop's edges: `10⁻⁶` of the body's size, at least +/// `10⁻⁵` mm — ten times [`planarity_tolerance`]. Truck meshes a face from its edges' own +/// polylines (shared with the neighbouring faces), so the mesh stays watertight whatever the +/// gap: measured on a five-sided rim, fills whose edges were up to `4.5·10⁻⁴` mm off the surface +/// meshed cleanly and stitched back into closed solids. Ten times the linear tolerance because a +/// fitted surface approaches curved edges like a cubic spline, so each further factor of ten +/// costs about twice the spans each way, and imported edges may already lie `10⁻⁵` mm or more off +/// their vertices; two orders below what meshing tolerates. Later operations must take the face +/// to contain its edges within this tolerance, not the linear one: the face's edges are the +/// body's exact edges, the surface between them is the fit. +pub(crate) fn face_tolerance(body: &Body) -> f64 { + (body.size() * 1e-6).max(1e-5) +} + +#[cfg(test)] +mod tests { + use mt::{BoundedCurve, InnerSpace, MetricSpace, ParametricCurve, ParametricSurface}; + use solidcraft_geom::{Loop2, Plane, Region2, Vec2, Vec3}; + use truck_modeling as mt; + + use crate::body::{Body, from_p3, p3}; + use crate::boolean_record::tests::{Reach, checked}; + use crate::history::{Change, EdgeOrigin}; + use crate::{BoolOp, boolean, copy_faces, cylinder, extrude, measure, patch_edges, stitch}; + + /// The cutting cylinder: axis along y at x = 0, z = 50, radius 40, so the top of a prism + /// within 20 of the axis and 20 high is cut away everywhere (to between 10 and 15.4 mm). + fn cut_height(x: f64) -> f64 { + 50.0 - (1600.0 - x * x).sqrt() + } + + /// A prism on the polygon `corners`, 20 high, its top cut by the cylinder; the cut solid, the + /// open cup without its top, and a point on each rim edge (the edges' midpoints in plan). + /// `inside` is a point inside the polygon. + fn cup(corners: &[Vec2], inside: Vec2) -> (Body, Body, Vec) { + let prism = extrude(&Plane::XY, &[Region2 { outer: Loop2::polygon(corners), holes: vec![] }], 0.0, 20.0).unwrap().pop().unwrap(); + let tool = cylinder(Vec3::new(0.0, -60.0, 50.0), Vec3::Y, 40.0, 120.0).unwrap(); + let solid = boolean(&prism, &tool, BoolOp::Cut).unwrap().unwrap(); + let n = corners.len(); + let mut keep = vec![inside.extend(0.0)]; + let mut rim = Vec::new(); + for k in 0..n { + let m = corners[k].lerp(corners[(k + 1) % n], 0.5); + keep.push(m.extend(5.0)); + rim.push(m.extend(cut_height(m.x))); + } + let open = copy_faces(&solid, &keep).unwrap(); + assert!(open.is_surface()); + (solid, open, rim) + } + + fn regular(n: usize, phase: f64) -> Vec { + (0..n).map(|k| Vec2::from_angle(phase + k as f64 * std::f64::consts::TAU / n as f64) * 20.0).collect() + } + + /// The distance from `p` to the surface: the nearest of a grid of samples, then Gauss–Newton + /// within the domain. + fn to_surface(s: &mt::Surface, p: mt::Point3) -> f64 { + let (Some((u0, u1)), Some((v0, v1))) = s.try_range_tuple() else { panic!("a bounded surface") }; + let n = 80; + let mut best = (f64::MAX, u0, v0); + for i in 0..=n { + for j in 0..=n { + let (u, v) = (u0 + (u1 - u0) * i as f64 / n as f64, v0 + (v1 - v0) * j as f64 / n as f64); + let d = s.subs(u, v).distance(p); + if d < best.0 { + best = (d, u, v); + } + } + } + let (_, mut u, mut v) = best; + for _ in 0..50 { + let (q, a, b) = (s.subs(u, v), s.uder(u, v), s.vder(u, v)); + let r = q - p; + let (aa, ab, bb) = (a.dot(a), a.dot(b), b.dot(b)); + let det = aa * bb - ab * ab; + if det.abs() < 1e-300 { + break; + } + let (ga, gb) = (a.dot(r), b.dot(r)); + u = (u - (bb * ga - ab * gb) / det).clamp(u0, u1); + v = (v - (aa * gb - ab * ga) / det).clamp(v0, v1); + } + s.subs(u, v).distance(p) + } + + /// The farthest point along the face's edges from its surface. + fn edges_off_surface(face: &mt::Face) -> f64 { + let s = face.surface(); + let mut worst = 0.0f64; + for e in face.boundaries().iter().flat_map(|w| w.edge_iter()) { + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + for k in 0..=24 { + worst = worst.max(to_surface(&s, c.subs(t0 + (t1 - t0) * k as f64 / 24.0))); + } + } + worst + } + + /// The fill's mesh: every normal finite, and no triangle folded over (all face the same way + /// as the mesh's total area vector, by a clear margin). + fn mesh_is_clean(fill: &Body) -> f64 { + let m = fill.tessellate(0.01).unwrap(); + assert!(!m.triangles.is_empty()); + assert!(m.normals.iter().all(|n| n.is_finite() && n.len() > 0.5), "a mesh normal is not finite"); + let tri = |t: &[u32; 3]| m.tri(t).unwrap(); + let total = m.triangles.iter().map(|t| { + let [a, b, c] = tri(t); + (b - a).cross(c - a) + }); + let total = total.fold(Vec3::ZERO, |s, v| s + v); + let n = total.normalized().unwrap(); + for t in &m.triangles { + let [a, b, c] = tri(t); + let v = (b - a).cross(c - a); + if v.len() > 1e-12 { + assert!(v.dot(n) > 0.0, "a folded triangle"); + } + } + measure(fill).unwrap().area + } + + /// Fill the rim, check the face, stitch it back: (the fill, the closed solid's volume). + fn fill_rim(open: &Body, rim: &[Vec3]) -> (Body, f64) { + let fill = patch_edges(open, rim).unwrap(); + assert!(fill.is_surface()); + assert_eq!(fill.face_count(), 1, "one face"); + let face = fill.solid.face_iter().next().unwrap().clone(); + assert!(!matches!(face.surface(), mt::Surface::Plane(_))); + // Its edges are the open body's own, by identity, and nothing else. + let own: Vec = open.solid.edge_iter().map(|e| e.id()).collect(); + assert!(fill.solid.edge_iter().all(|e| own.contains(&e.id())), "an edge that is not the body's"); + let h = checked(&fill, &[open], 0.05, Reach { made: 0.0, moved: 0.0 }, &|_| true); + assert!(h.faces.iter().all(|f| f.change == Change::Generated && f.origins.is_empty())); + assert!(h.edges.iter().all(|e| matches!(e.origins.as_slice(), [EdgeOrigin::Edge { operand: 0, .. }]))); + let off = edges_off_surface(&face); + let tol = super::face_tolerance(open); + assert!(off <= tol, "an edge is {off} off the face (tolerance {tol})"); + mesh_is_clean(&fill); + let closed = stitch(&[open, &fill], 0.0).unwrap(); + assert!(!closed.is_surface(), "the fill closes the cup"); + (fill, measure(&closed).unwrap().volume) + } + + /// The fill's height over the point `(x, y)`: Newton's method on the surface's `x`, `y` from + /// the nearest of a grid of samples. + fn height_over(s: &mt::Surface, x: f64, y: f64, start: &mut Option<(f64, f64)>) -> Option { + let (Some((u0, u1)), Some((v0, v1))) = s.try_range_tuple() else { return None }; + let (mut u, mut v) = start.unwrap_or(((u0 + u1) / 2.0, (v0 + v1) / 2.0)); + for _ in 0..30 { + let (q, a, b) = (s.subs(u, v), s.uder(u, v), s.vder(u, v)); + let (rx, ry) = (q.x - x, q.y - y); + let det = a.x * b.y - a.y * b.x; + if det.abs() < 1e-300 { + return None; + } + u = (u - (b.y * rx - b.x * ry) / det).clamp(u0, u1); + v = (v - (a.x * ry - a.y * rx) / det).clamp(v0, v1); + } + let q = s.subs(u, v); + *start = Some((u, v)); + ((q.x - x).hypot(q.y - y) < 1e-9).then_some(q.z) + } + + /// The volume the cut prism gains with the fill for its top: the fill's height above the + /// cylinder integrated over the polygon (midpoint rule on a 300 × 300 grid; the integrand + /// vanishes on the rim, so cells across it add little error). + fn gained(fill: &Body, corners: &[Vec2]) -> f64 { + let s = fill.solid.face_iter().next().unwrap().surface(); + let (lo, hi) = corners.iter().fold((Vec2::new(f64::MAX, f64::MAX), Vec2::new(f64::MIN, f64::MIN)), |(lo, hi), p| { + (Vec2::new(lo.x.min(p.x), lo.y.min(p.y)), Vec2::new(hi.x.max(p.x), hi.y.max(p.y))) + }); + let n = 300; + let (dx, dy) = ((hi.x - lo.x) / n as f64, (hi.y - lo.y) / n as f64); + let inside = |x: f64, y: f64| { + let mut w = false; + for k in 0..corners.len() { + let (a, b) = (corners[k], corners[(k + 1) % corners.len()]); + if (a.y > y) != (b.y > y) && x < a.x + (y - a.y) / (b.y - a.y) * (b.x - a.x) { + w = !w; + } + } + w + }; + let mut sum = 0.0; + let mut start = None; + for i in 0..n { + let x = lo.x + (i as f64 + 0.5) * dx; + for j in 0..n { + let y = lo.y + (j as f64 + 0.5) * dy; + if inside(x, y) { + let h = height_over(&s, x, y, &mut start).expect("the fill covers the polygon"); + sum += (h - cut_height(x)) * dx * dy; + } + } + } + sum + } + + /// The closed solid's volume against the cut prism's plus what the fill adds over the + /// cylinder (integrated independently, [`gained`]); they agree to the integration's error + /// (0.1% of what the fill adds) plus the error of measuring both solids on meshes + /// (`5·10⁻⁵` of the volume). + fn check_volume(name: &str, solid: &Body, fill: &Body, corners: &[Vec2], got: f64) { + let prism = measure(solid).unwrap().volume; + let extra = gained(fill, corners); + let want = prism + extra; + let bound = 1e-3 * extra.abs() + 5e-5 * want; + eprintln!("{name}: closed volume {got:.4}, cut prism {prism:.4} + fill above the cylinder {extra:.4} = {want:.4} (bound ±{bound:.4})"); + assert!((got - want).abs() <= bound, "{name}: {got} vs {want}"); + } + + #[test] + fn regular_rims_of_three_to_eight_sides_fill_with_one_face_and_close_the_cup() { + for n in 3..=8 { + let corners = regular(n, 0.3); + let (solid, open, rim) = cup(&corners, Vec2::ZERO); + let (fill, got) = fill_rim(&open, &rim); + assert_eq!(fill.edge_count(), n, "the rim's own edges"); + check_volume(&format!("{n} sides"), &solid, &fill, &corners, got); + } + } + + #[test] + fn concave_rims_fill_with_one_face() { + let v = |x: f64, y: f64| Vec2::new(x, y); + let shapes = [ + ("dart", vec![v(-10.0, 0.0), v(15.0, -12.0), v(3.0, 0.0), v(15.0, 12.0)], v(0.0, 0.0)), + ("arrow", vec![v(-12.0, 0.0), v(0.0, -10.0), v(0.0, -4.0), v(14.0, 0.0), v(0.0, 4.0)], v(0.0, 0.0)), + ( + "plus", + vec![ + v(-4.0, -14.0), + v(4.0, -14.0), + v(4.0, -4.0), + v(14.0, -4.0), + v(14.0, 4.0), + v(4.0, 4.0), + v(4.0, 14.0), + v(-4.0, 14.0), + v(-4.0, 4.0), + v(-14.0, 4.0), + v(-14.0, -4.0), + v(-4.0, -4.0), + ], + v(0.0, 0.0), + ), + ]; + for (name, corners, inside) in shapes { + let (solid, open, rim) = cup(&corners, inside); + let (fill, got) = fill_rim(&open, &rim); + assert_eq!(fill.edge_count(), corners.len(), "{name}: the rim's own edges"); + check_volume(name, &solid, &fill, &corners, got); + } + } + + #[test] + fn a_flat_corner_between_two_edges_is_kept() { + // A rim of five straight edges, two of them in line (a 180° corner between them), the + // four true corners at different heights; walls hang below it. + let p = |x: f64, y: f64, z: f64| mt::Vertex::new(mt::Point3::new(x, y, z)); + let v = [p(-14.0, -14.0, 12.0), p(0.0, -14.0, 12.0), p(14.0, -14.0, 12.0), p(14.0, 14.0, 15.0), p(-14.0, 14.0, 10.0)]; + let line = |a: &mt::Vertex, b: &mt::Vertex| mt::Edge::new(a, b, mt::Curve::Line(mt::Line(a.point(), b.point()))); + let wire: mt::Wire = (0..5).map(|k| line(&v[k], &v[(k + 1) % 5])).collect::>().into(); + let low = mt::builder::translated(&wire, mt::Vector3::new(0.0, 0.0, -10.0)); + let walls = mt::builder::try_wire_homotopy(&wire, &low).unwrap(); + let body = crate::surfaces::surface_body(walls.face_iter().cloned().collect()); + let rim: Vec = (0..5).map(|k| from_p3(v[k].point()).lerp(from_p3(v[(k + 1) % 5].point()), 0.5)).collect(); + let fill = patch_edges(&body, &rim).unwrap(); + assert_eq!((fill.face_count(), fill.edge_count()), (1, 5), "one face on all five edges"); + let face = fill.solid.face_iter().next().unwrap().clone(); + assert!(edges_off_surface(&face) <= super::face_tolerance(&body)); + let area = mesh_is_clean(&fill); + eprintln!("flat corner: one face of area {area}"); + assert!(area > 28.0 * 28.0 * 0.99 && area < 28.0 * 28.0 * 1.05, "{area}"); + } + + #[test] + fn the_fill_does_not_depend_on_pick_order_or_where_the_body_is() { + let (_, open, rim) = cup(®ular(5, 0.3), Vec2::ZERO); + let (a, _) = fill_rim(&open, &rim); + let mut shuffled = rim.clone(); + shuffled.rotate_left(3); + shuffled.swap(0, 1); + let (b, _) = fill_rim(&open, &shuffled); + let (sa, sb) = (a.solid.face_iter().next().unwrap().surface(), b.solid.face_iter().next().unwrap().surface()); + let (Some((u0, u1)), Some((v0, v1))) = sa.try_range_tuple() else { panic!() }; + let mut gap = 0.0f64; + for (x, y) in [(0.2, 0.3), (0.5, 0.5), (0.7, 0.9)] { + let (u, v) = (u0 + (u1 - u0) * x, v0 + (v1 - v0) * y); + gap = gap.max(sa.subs(u, v).distance(sb.subs(u, v))); + } + assert!(gap < 1e-7, "pick order moved the surface by {gap}"); + let (area_a, area_b) = (mesh_is_clean(&a), mesh_is_clean(&b)); + assert!((area_a - area_b).abs() < area_a * 1e-9); + // Moved and turned: the same area. + let moved = crate::transform(&open, Vec3::new(130.0, -70.0, 40.0), Vec3::ZERO, Vec3::new(1.0, 2.0, 3.0), 0.7).unwrap(); + let ax = Vec3::new(1.0, 2.0, 3.0).normalized().unwrap(); + let turn = |p: Vec3| { + // Rodrigues' rotation by 0.7 about `ax`, then the shift. + let (c, s) = (0.7f64.cos(), 0.7f64.sin()); + p * c + ax.cross(p) * s + ax * (ax.dot(p) * (1.0 - c)) + Vec3::new(130.0, -70.0, 40.0) + }; + let moved_rim: Vec = rim.iter().map(|p| turn(*p)).collect(); + let (c, _) = fill_rim(&moved, &moved_rim); + let area_c = mesh_is_clean(&c); + eprintln!("pick order: surfaces within {gap:e}; areas {area_a} / {area_b} / moved {area_c}"); + assert!((area_a - area_c).abs() < area_a * 1e-6, "{area_a} vs moved {area_c}"); + } + + #[test] + fn a_saddle_rim_of_two_intersection_curves_fills() { + // A round post with its top cut by the sideways cylinder: the rim is two intersection + // curves; the fill is one face on them. + let post = cylinder(Vec3::ZERO, Vec3::Z, 10.0, 20.0).unwrap(); + let tool = cylinder(Vec3::new(0.0, -60.0, 50.0), Vec3::Y, 40.0, 120.0).unwrap(); + let solid = boolean(&post, &tool, BoolOp::Cut).unwrap().unwrap(); + let sides: Vec = (0..8).map(|k| Vec2::from_angle(k as f64 * 0.8).extend(0.0) * 10.0 + Vec3::new(0.0, 0.0, 5.0)).collect(); + let open = copy_faces(&solid, &[&[Vec3::ZERO][..], &sides].concat()).unwrap(); + let rim: Vec = [-0.98, 2.18] + .map(|a: f64| { + let (x, y) = (10.0 * a.cos(), 10.0 * a.sin()); + Vec3::new(x, y, cut_height(x)) + }) + .to_vec(); + let (fill, got) = fill_rim(&open, &rim); + let circle: Vec = (0..720).map(|k| Vec2::from_angle(k as f64 * std::f64::consts::TAU / 720.0) * 10.0).collect(); + check_volume("saddle rim", &solid, &fill, &circle, got); + } + + #[test] + fn imported_edges_whose_vertices_are_off_their_curves_fill() { + // The cup's vertices moved 2·10⁻⁵ mm off their curves, as a file might give them. + let corners = regular(5, 0.3); + let (solid, open, rim) = cup(&corners, Vec2::ZERO); + let shell = open.solid.boundaries()[0].mapped( + |p: &mt::Point3| *p + mt::Vector3::new(2e-5, 0.0, 0.0), + |c: &mt::Curve| c.clone(), + |s: &mt::Surface| s.clone(), + ); + let imported = crate::surfaces::surface_body(shell.face_iter().cloned().collect()); + let (fill, got) = fill_rim(&imported, &rim); + check_volume("imported", &solid, &fill, &corners, got); + } + + #[test] + fn a_spline_edge_bulging_out_of_plane_is_never_given_a_flat_face() { + // A flat rectangle of edges but for one cubic edge with a control point raised 1.5 mm: no + // plane can be certified. A one-span bump is followed by the fitted face; a narrow one + // (64 spans, about 1 mm wide, its third derivative jumping at every knot) may be beyond + // a 32 × 32 grid, and then it is refused. + for (spans, raised, must_fill) in [(1, 1, true), (64, 30, false)] { + let count = spans + 3; + let pts: Vec = + (0..count).map(|i| mt::Point3::new(20.0 * i as f64 / (count - 1) as f64, 0.0, if i == raised { 1.5 } else { 0.0 })).collect(); + let spline = mt::BSplineCurve::new(mt::KnotVec::uniform_knot(3, spans), pts); + let v: Vec = + [(0.0, 0.0), (20.0, 0.0), (20.0, 10.0), (0.0, 10.0)].iter().map(|&(x, y)| mt::Vertex::new(mt::Point3::new(x, y, 0.0))).collect(); + let line = |a: &mt::Vertex, b: &mt::Vertex| mt::Edge::new(a, b, mt::Curve::Line(mt::Line(a.point(), b.point()))); + let bottom = mt::Edge::new(&v[0], &v[1], mt::Curve::BSplineCurve(spline.clone())); + let wire: mt::Wire = vec![bottom, line(&v[1], &v[2]), line(&v[2], &v[3]), line(&v[3], &v[0])].into(); + let top = mt::builder::translated(&wire, mt::Vector3::new(0.0, 0.0, -5.0)); + let walls = mt::builder::try_wire_homotopy(&wire, &top).unwrap(); + let body = crate::surfaces::surface_body(walls.face_iter().cloned().collect()); + let peak = from_p3(spline.subs((raised as f64 + 1.0) / (count + 1) as f64)); + let rim = [from_p3(spline.subs(0.5)), Vec3::new(20.0, 5.0, 0.0), Vec3::new(10.0, 10.0, 0.0), Vec3::new(0.0, 5.0, 0.0)]; + match patch_edges(&body, &rim) { + Ok(fill) => { + let face = fill.solid.face_iter().next().unwrap().clone(); + assert!(!matches!(face.surface(), mt::Surface::Plane(_)), "a plane would miss the bulge"); + let (off, tol) = (edges_off_surface(&face), super::face_tolerance(&body)); + eprintln!("{spans}-span bump: filled, edges within {off:e} of the face"); + assert!(off <= tol && to_surface(&face.surface(), p3(peak)) <= tol); + mesh_is_clean(&fill); + } + Err(e) => { + eprintln!("{spans}-span bump: {e}"); + assert!(!must_fill && e.to_string().contains("not supported yet"), "{e}"); + } + } + } + } + + /// A point on each edge that bounds one face of `body` (its free edges). + fn free_edge_points(body: &Body) -> Vec { + let mut uses: std::collections::HashMap = std::collections::HashMap::new(); + for f in body.solid.face_iter() { + for e in f.boundaries().iter().flat_map(|w| w.edge_iter()) { + uses.entry(e.id()).or_insert((0, e.clone())).0 += 1; + } + } + uses.values() + .filter(|(n, _)| *n == 1) + .map(|(_, e)| { + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + from_p3(c.subs(0.5 * (t0 + t1))) + }) + .collect() + } + + #[test] + fn a_loop_that_does_not_bound_a_height_field_is_refused() { + // Half a round post's side: two straight edges up and two half circles, whose plane + // stands on the post's axis; the half circles run along its normal at their ends. + let post = cylinder(Vec3::ZERO, Vec3::Z, 10.0, 20.0).unwrap(); + let halves = crate::split_by_plane(&post, &Plane::from_normal(Vec3::ZERO, Vec3::X).unwrap()).unwrap(); + let half = halves.into_iter().find(|b| measure(b).unwrap().bbox.min.x > -1e-6).unwrap(); + let at: Vec = + (-4..=4).map(|k| Vec2::from_angle(k as f64 * 0.35).extend(10.0) * 1.0).map(|p| Vec3::new(p.x * 10.0, p.y * 10.0, 10.0)).collect(); + let side = copy_faces(&half, &at).unwrap(); + let rim = free_edge_points(&side); + let err = patch_edges(&side, &rim).unwrap_err().to_string(); + eprintln!("half a post's side: {err}"); + assert!(err.contains("not supported yet"), "{err}"); + } + + #[test] + fn planar_loops_of_intersection_curves_are_one_planar_face() { + // A round post cut by planes tilted every way, gently to steeply: each rim is an ellipse + // of intersection curves (split where the post's side has seams) lying on the cutting + // plane, so the patch is one planar face, and the post closes again. + for tilt in [Vec3::new(0.0, -0.4, 1.0), Vec3::new(0.7, -0.9, 1.0), Vec3::new(1.3, 0.2, 1.0)] { + let post = cylinder(Vec3::ZERO, Vec3::Z, 10.0, 30.0).unwrap(); + let n = tilt.normalized().unwrap(); + let plane = Plane::from_normal(Vec3::new(0.0, 0.0, 15.0), n).unwrap(); + let pieces = crate::split_by_plane(&post, &plane).unwrap(); + let solid = pieces.into_iter().find(|b| measure(b).unwrap().bbox.min.z < 0.5).unwrap(); + let height = |x: f64, y: f64| 15.0 - (tilt.x * x + tilt.y * y) / tilt.z; + let sides: Vec = (0..8).map(|k| Vec2::from_angle(k as f64 * 0.8).extend(0.0) * 10.0 + Vec3::new(0.0, 0.0, 1.0)).collect(); + let open = copy_faces(&solid, &[&[Vec3::ZERO][..], &sides].concat()).unwrap(); + let rim: Vec = (0..12) + .map(|k| { + let a = k as f64 * std::f64::consts::TAU / 12.0 + 0.2; + let (x, y) = (10.0 * a.cos(), 10.0 * a.sin()); + Vec3::new(x, y, height(x, y)) + }) + .collect(); + let lid = patch_edges(&open, &rim).unwrap(); + assert_eq!(lid.face_count(), 1); + let edges = lid.edge_count(); + assert!(matches!(lid.solid.face_iter().next().unwrap().surface(), mt::Surface::Plane(_)), "tilt {tilt:?}"); + let closed = stitch(&[&open, &lid], 0.0).unwrap(); + assert!(!closed.is_surface()); + let (want, got) = (measure(&solid).unwrap().volume, measure(&closed).unwrap().volume); + eprintln!("rim cut at {tilt:?}: {edges} intersection curves, planar face, volume {got} against {want}"); + assert!((got - want).abs() < want * 1e-6, "{got} vs {want}"); + } + } + + #[test] + fn planar_loops_are_still_one_planar_face() { + let b = crate::box_solid(Vec3::ZERO, Vec3::new(10.0, 10.0, 10.0)).unwrap(); + let open = copy_faces( + &b, + &[Vec3::new(5.0, 5.0, 0.0), Vec3::new(5.0, 0.0, 5.0), Vec3::new(10.0, 5.0, 5.0), Vec3::new(5.0, 10.0, 5.0), Vec3::new(0.0, 5.0, 5.0)], + ) + .unwrap(); + let rim = [Vec3::new(5.0, 0.0, 10.0), Vec3::new(10.0, 5.0, 10.0), Vec3::new(5.0, 10.0, 10.0), Vec3::new(0.0, 5.0, 10.0)]; + let lid = patch_edges(&open, &rim).unwrap(); + assert_eq!(lid.face_count(), 1); + assert!(matches!(lid.solid.face_iter().next().unwrap().surface(), mt::Surface::Plane(_))); + let closed = stitch(&[&open, &lid], 0.0).unwrap(); + assert!((measure(&closed).unwrap().volume - 1000.0).abs() < 1e-9); + } + + #[test] + fn open_partial_or_hostile_picks_are_refused() { + let (_, open, rim) = cup(®ular(5, 0.3), Vec2::ZERO); + assert!(patch_edges(&open, &rim[..4]).is_err(), "a loop with a missing edge"); + assert!(patch_edges(&open, &[]).is_err()); + assert!(patch_edges(&open, &[Vec3::new(f64::NAN, 0.0, 0.0)]).is_err()); + assert!(patch_edges(&open, &[Vec3::new(1e300, 0.0, 0.0)]).is_err()); + assert!(patch_edges(&open, &vec![rim[0]; 20_000]).is_err()); + } +} diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index ec231d65..7defd097 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -26,6 +26,7 @@ mod delete_face; mod emboss; #[cfg(test)] mod face_source_tests; +mod fill; mod freeform; pub mod fuzz; mod heal; @@ -45,6 +46,7 @@ mod modify_history_tests; mod numbered; mod offset; mod ops; +mod plane_fit; mod polybool; mod polyhedron; mod prism; @@ -54,6 +56,7 @@ pub mod provenance; mod record; mod sew; mod smooth_edges; +mod spline_conv; mod splitface; mod step; mod step_in; diff --git a/crates/kernel/src/plane_fit.rs b/crates/kernel/src/plane_fit.rs new file mode 100644 index 00000000..bfb118bd --- /dev/null +++ b/crates/kernel/src/plane_fit.rs @@ -0,0 +1,133 @@ +//! Planes fitted to points around a loop (Newell's method), for planar faces of polygons read +//! from files and for patches filling a planar loop of edges, and points sampled along edges +//! (shared with the fill of a loop that is not planar). + +use mt::{BoundedCurve, EuclideanSpace, InnerSpace, MetricSpace, ParametricCurve}; +use solidcraft_geom::Vec3; +use truck_modeling as mt; + +use crate::body::{from_p3, p3}; +use crate::{KernelError, Result}; + +/// The plane through points round a loop: its normal by Newell's method (turning with the +/// loop), through their centroid; `None` for fewer than three points or no area. +pub(crate) fn plane_through(pts: &[mt::Point3]) -> Option { + let n = pts.len(); + if n < 3 { + return None; + } + let mut nrm = mt::Vector3::new(0.0, 0.0, 0.0); + let mut c = mt::Vector3::new(0.0, 0.0, 0.0); + for i in 0..n { + let (a, b) = (pts.get(i)?, pts.get((i + 1) % n)?); + nrm.x += (a.y - b.y) * (a.z + b.z); + nrm.y += (a.z - b.z) * (a.x + b.x); + nrm.z += (a.x - b.x) * (a.y + b.y); + c += a.to_vec(); + } + let len = nrm.magnitude(); + if !(len > 1e-12) { + return None; + } + let z = nrm / len; + let o = mt::Point3::from_vec(c / n as f64); + let a = if z.x.abs() < 0.9 { mt::Vector3::unit_x() } else { mt::Vector3::unit_y() }; + let x = (a - z * a.dot(z)).normalize(); + let y = z.cross(x); + Some(mt::Plane::new(o, o + x, o + y)) +} + +/// The point of `c` at `t` (`None` where an intersection curve cannot find it). +fn curve_point(c: &mt::Curve, t: f64) -> Option { + let p = match c { + mt::Curve::IntersectionCurve(ic) => ic.search_triple(t)?.0, + other => other.subs(t), + }; + let p = from_p3(p); + p.is_finite().then_some(p) +} + +/// Points along an edge, in its own direction, at `s ∈ [0, 1]` of its parameter range. An +/// intersection curve's point is searched from the surface parameters of the point found before +/// when that one is near (an eighth of the range), and kept only when it lies near the curve's +/// leader (a seeded search can slide to another branch of the surfaces' intersection); else the +/// unseeded search decides. `None` where no point is found. +pub(crate) struct EdgeSampler { + curve: mt::Curve, + range: (f64, f64), + hint: Option<(f64, mt::Point2, mt::Point2)>, +} + +impl EdgeSampler { + pub(crate) fn new(edge: &mt::Edge) -> Option { + let curve = edge.oriented_curve(); + let range = curve.range_tuple(); + (range.0.is_finite() && range.1.is_finite() && range.0 < range.1).then_some(EdgeSampler { curve, range, hint: None }) + } + + pub(crate) fn at(&mut self, s: f64) -> Option { + let (t0, t1) = self.range; + let t = if s >= 1.0 { t1 } else { t0 + (t1 - t0) * s }; + let mt::Curve::IntersectionCurve(ic) = &self.curve else { return curve_point(&self.curve, t) }; + let leader = ic.leader(); + let (p, d) = (leader.subs(t), leader.der(t)); + let near = |q: mt::Point3| q.distance(p) <= 4.0 * ic.tolerance() + 1e-6; + let seeded = self + .hint + .filter(|(tp, _, _)| (t - tp).abs() <= (t1 - t0) / 8.0) + .and_then(|(_, a, b)| mt::double_projection(ic.surface0(), Some((a.x, a.y)), ic.surface1(), Some((b.x, b.y)), p, d, 100)) + .filter(|f| near(f.0)); + let found = seeded.or_else(|| mt::double_projection(ic.surface0(), None, ic.surface1(), None, p, d, 100))?; + self.hint = Some((t, found.1, found.2)); + let q = from_p3(found.0); + q.is_finite().then_some(q) + } +} + +/// The planar face bounded by the closed loop `chain` (each edge oriented along it) when its +/// edges are certainly in one plane, `None` otherwise (also when an edge cannot be sampled). +/// The plane is fitted to points along the edges ([`plane_through`], the normal turning with the +/// loop); the edges are certified within `tol` of it conservatively: a line by its ends, a +/// B-spline or NURBS curve with positive weights by its control points (the curve lies in their +/// convex hull), an intersection curve only when one of its two surfaces is a plane that the +/// loop's points lie on and the fitted plane agrees with over the loop (the curve lies in that +/// plane). Samples alone are never trusted: a spline can bulge out of the plane between them. +pub(crate) fn planar_face(chain: &[mt::Edge], tol: f64) -> Result> { + const SAMPLES: usize = 16; + let mut pts = Vec::with_capacity(chain.len() * SAMPLES); + for e in chain { + let Some(mut sampler) = EdgeSampler::new(e) else { return Ok(None) }; + for k in 0..SAMPLES { + let Some(p) = sampler.at(k as f64 / SAMPLES as f64) else { return Ok(None) }; + pts.push(p3(p)); + } + } + let Some(plane) = plane_through(&pts) else { return Ok(None) }; + let (o, n) = (plane.origin(), plane.normal()); + let off = |p: mt::Point3, o: mt::Point3, n: mt::Vector3| (p - o).dot(n).abs(); + let on_plane = |p: mt::Point3| off(p, o, n) <= tol; + if !pts.iter().all(|p| on_plane(*p)) { + return Ok(None); + } + // A support plane `q` of an intersection curve holds the loop when the loop's points lie on + // it and their feet on the fitted plane do too. + let holds = |q: &mt::Plane| { + let (qo, qn) = (q.origin(), q.normal()); + pts.iter().all(|p| off(*p, qo, qn) <= tol && off(*p - n * (*p - o).dot(n), qo, qn) <= tol) + }; + for e in chain { + let certified = match e.curve() { + mt::Curve::Line(l) => on_plane(l.0) && on_plane(l.1), + mt::Curve::BSplineCurve(b) => b.control_points().iter().all(|p| on_plane(*p)), + mt::Curve::NurbsCurve(c) => { + c.non_rationalized().control_points().iter().all(|h| h.w > 0.0 && on_plane(mt::Point3::new(h.x / h.w, h.y / h.w, h.z / h.w))) + } + mt::Curve::IntersectionCurve(ic) => [ic.surface0(), ic.surface1()].into_iter().any(|s| matches!(s, mt::Surface::Plane(q) if holds(q))), + }; + if !certified { + return Ok(None); + } + } + let wire: mt::Wire = chain.to_vec().into(); + Ok(Some(mt::Face::try_new(vec![wire], mt::Surface::Plane(plane)).map_err(|e| KernelError::Failed(format!("patch: {e}")))?)) +} diff --git a/crates/kernel/src/spline_conv.rs b/crates/kernel/src/spline_conv.rs new file mode 100644 index 00000000..5a681220 --- /dev/null +++ b/crates/kernel/src/spline_conv.rs @@ -0,0 +1,36 @@ +//! geom's B-spline curves and surfaces as truck's NURBS geometry, weights kept (homogeneous +//! control points carried over as they are). + +use solidcraft_geom::bspline::{HPoint, NurbsSurface}; +use truck_modeling as mt; + +use crate::{KernelError, Result}; + +fn fail(m: impl std::fmt::Display) -> KernelError { + KernelError::Failed(format!("spline geometry: {m}")) +} + +/// A homogeneous control point as truck's. +pub(crate) fn h4(h: HPoint) -> mt::Vector4 { + mt::Vector4::new(h.x, h.y, h.z, h.w) +} + +/// A NURBS curve of knots and homogeneous control points. +pub(crate) fn nurbs_curve(knots: mt::KnotVec, ctrl: &[HPoint]) -> Result { + let c = mt::BSplineCurve::try_new(knots, ctrl.iter().copied().map(h4).collect()).map_err(fail)?; + Ok(mt::Curve::NurbsCurve(mt::NurbsCurve::new(c))) +} + +/// A NURBS surface of knots and homogeneous control points (`rows[i][j]`: i along u). +pub(crate) fn nurbs_surface(knots: (mt::KnotVec, mt::KnotVec), rows: Vec>) -> Result { + let s = mt::BSplineSurface::try_new(knots, rows).map_err(fail)?; + Ok(mt::Surface::NurbsSurface(mt::NurbsSurface::new(s))) +} + +/// geom's tensor surface as truck's. +pub(crate) fn surface(s: &NurbsSurface) -> Result { + let [nu, nv] = s.count(); + let [ku, kv] = s.knots(); + let rows: Vec> = (0..nu).map(|i| (0..nv).map(|j| h4(s.point(i, j))).collect()).collect(); + nurbs_surface((mt::KnotVec::from(ku.clone()), mt::KnotVec::from(kv.clone())), rows) +} diff --git a/crates/kernel/src/step_in/brep.rs b/crates/kernel/src/step_in/brep.rs index cf8e4b30..33f110eb 100644 --- a/crates/kernel/src/step_in/brep.rs +++ b/crates/kernel/src/step_in/brep.rs @@ -4,7 +4,7 @@ use std::collections::HashMap; -use mt::{BoundedCurve, EuclideanSpace, InnerSpace, MetricSpace, ParametricCurve}; +use mt::{BoundedCurve, InnerSpace, MetricSpace, ParametricCurve}; use truck_modeling as mt; use super::Ctx; @@ -244,7 +244,7 @@ impl Builder { }) .collect(); let surface = if e.name() == "FACE" { - plane_through(&extent).ok_or("polygonal face without a plane")? + crate::plane_fit::plane_through(&extent).map(mt::Surface::Plane).ok_or("polygonal face without a plane")? } else { let sid = rid(p.get(2), "face geometry")?; let sense = p.get(3).and_then(Param::as_bool).unwrap_or(true); @@ -344,30 +344,6 @@ impl Builder { } } -/// Plane through a polygon (Newell's normal). -fn plane_through(pts: &[mt::Point3]) -> Option { - let n = pts.len(); - let mut nrm = mt::Vector3::new(0.0, 0.0, 0.0); - let mut c = mt::Vector3::new(0.0, 0.0, 0.0); - for i in 0..n { - let (a, b) = (pts.get(i)?, pts.get((i + 1) % n)?); - nrm.x += (a.y - b.y) * (a.z + b.z); - nrm.y += (a.z - b.z) * (a.x + b.x); - nrm.z += (a.x - b.x) * (a.y + b.y); - c += a.to_homogeneous().truncate(); - } - let len = nrm.magnitude(); - if !(len > 1e-12) || n == 0 { - return None; - } - let z = nrm / len; - let o = mt::Point3::from_vec(c / n as f64); - let a = if z.x.abs() < 0.9 { mt::Vector3::unit_x() } else { mt::Vector3::unit_y() }; - let x = (a - z * a.dot(z)).normalize(); - let y = z.cross(x); - Some(mt::Surface::Plane(mt::Plane::new(o, o + x, o + y))) -} - /// Read a solid-like representation item: MANIFOLD_SOLID_BREP, BREP_WITH_VOIDS, FACETED_BREP, /// or a SHELL_BASED_SURFACE_MODEL (closed shells give solids, open ones open bodies). pub(crate) fn read_item(cx: &Ctx, id: u64) -> R { diff --git a/crates/kernel/src/surfaces.rs b/crates/kernel/src/surfaces.rs index b1392482..366fc14f 100644 --- a/crates/kernel/src/surfaces.rs +++ b/crates/kernel/src/surfaces.rs @@ -1,5 +1,6 @@ -//! Surface bodies (open shells) and the operations that turn them into solids: Patch (a planar -//! face filling a closed loop of edges or a sketch region), Stitch (surfaces sharing edges sewn +//! Surface bodies (open shells) and the operations that turn them into solids: Patch (a face +//! filling a closed loop of edges — planar, or B-spline patches for a loop that is not — or a +//! planar sketch region), Stitch (surfaces sharing edges sewn //! into one body, a solid when they close), Trim (a surface cut by a tool, one side kept), //! Thicken (a surface offset along its normals into a solid) and copying faces out of a body. @@ -76,13 +77,31 @@ pub(crate) fn edge_near(body: &Body, p: Vec3) -> Option { best.filter(|(d, _)| *d <= body.size() * 0.05 + 1e-6).map(|(_, e)| e) } -/// A planar surface filling a closed loop of a body's edges (each picked by a point on it). Its -/// history is over `[body]`: the face is new, its edges the body's. +/// A surface filling a closed loop of a body's edges (each picked by a point on it): one face on +/// the body's own edges, planar when the loop is certainly planar ([`crate::plane_fit`]), else a +/// B-spline height field over the loop's plane ([`crate::fill`]; loops that do not bound one are +/// refused). Its history is over `[body]`: the face is new, its edges the body's own. pub fn patch_edges(body: &Body, at: &[Vec3]) -> Result { body.require_brep("patch")?; if at.is_empty() || at.len() > 10_000 { return Err(KernelError::Invalid("pick the edges of a closed loop".into())); } + guard("patch", || { + let chain = picked_loop(body, at)?; + let Some(face) = crate::plane_fit::planar_face(&chain, crate::fill::planarity_tolerance(body))? else { + // Not certainly planar: a fitted surface. + return crate::fill::fill_loop(body, &chain); + }; + // A new face on the body's own edges. + let mut rec = Rebuild::default(); + rec.made(&face, Vec::new()); + Ok(rebuilt_body(surface_body(vec![face]), &[&body.solid], &rec)) + }) +} + +/// The edges of `body` picked by points on them, chained into one closed loop (each turned to +/// run along it). Call it under [`guard`]: sampling an intersection curve can fail in truck. +pub(crate) fn picked_loop(body: &Body, at: &[Vec3]) -> Result> { let mut edges: Vec = Vec::new(); for p in at { let e = edge_near(body, *p).ok_or_else(|| KernelError::Invalid(format!("no edge at {p:?}")))?; @@ -90,34 +109,25 @@ pub fn patch_edges(body: &Body, at: &[Vec3]) -> Result { edges.push(e); } } - guard("patch", || { - // Chain the edges into one loop (turning them round as needed). - let mut pool = edges; - let first = pool.pop().ok_or_else(|| KernelError::Invalid("no edges".into()))?; - let start = first.front().id(); - let mut at_v = first.back().clone(); - let mut chain = vec![first]; - while at_v.id() != start { - let k = pool - .iter() - .position(|e| e.front().id() == at_v.id() || e.back().id() == at_v.id()) - .ok_or_else(|| KernelError::Invalid("the edges do not form a closed loop".into()))?; - let e = pool.remove(k); - let e = if e.front().id() == at_v.id() { e } else { e.inverse() }; - at_v = e.back().clone(); - chain.push(e); - } - if !pool.is_empty() { - return Err(KernelError::Invalid("the edges form more than one loop".into())); - } - let wire: mt::Wire = chain.into(); - let face = mt::builder::try_attach_plane(&[wire]) - .map_err(|_| KernelError::Invalid("not supported yet: patching a loop that is not planar".into()))?; - // A new face on the body's own edges. - let mut rec = Rebuild::default(); - rec.made(&face, Vec::new()); - Ok(rebuilt_body(surface_body(vec![face]), &[&body.solid], &rec)) - }) + let mut pool = edges; + let first = pool.pop().ok_or_else(|| KernelError::Invalid("no edges".into()))?; + let start = first.front().id(); + let mut at_v = first.back().clone(); + let mut chain = vec![first]; + while at_v.id() != start { + let k = pool + .iter() + .position(|e| e.front().id() == at_v.id() || e.back().id() == at_v.id()) + .ok_or_else(|| KernelError::Invalid("the edges do not form a closed loop".into()))?; + let e = pool.remove(k); + let e = if e.front().id() == at_v.id() { e } else { e.inverse() }; + at_v = e.back().clone(); + chain.push(e); + } + if !pool.is_empty() { + return Err(KernelError::Invalid("the edges form more than one loop".into())); + } + Ok(chain) } /// Faces of a body copied out as a surface body (faces picked by points on them). Its history is diff --git a/crates/kernel/src/tspline.rs b/crates/kernel/src/tspline.rs index 2a6e6a9b..aa864399 100644 --- a/crates/kernel/src/tspline.rs +++ b/crates/kernel/src/tspline.rs @@ -11,6 +11,7 @@ use solidcraft_geom::tspline::TSpline; use truck_modeling as mt; use crate::body::Body; +use crate::spline_conv::{h4, nurbs_curve, nurbs_surface}; use crate::{KernelError, Result, guard}; /// A T-spline as a surface body (an open shell of NURBS faces whose normals follow @@ -34,25 +35,11 @@ fn fail(m: impl std::fmt::Display) -> KernelError { KernelError::Failed(format!("t-spline sheet: {m}")) } -fn h4(h: HPoint) -> mt::Vector4 { - mt::Vector4::new(h.x, h.y, h.z, h.w) -} - fn vertex(h: HPoint) -> Result { let p = h.project().ok_or_else(|| fail("a corner has no position"))?; Ok(mt::builder::vertex(crate::body::p3(p))) } -fn nurbs_curve(knots: mt::KnotVec, ctrl: &[HPoint]) -> Result { - let c = mt::BSplineCurve::try_new(knots, ctrl.iter().copied().map(h4).collect()).map_err(fail)?; - Ok(mt::Curve::NurbsCurve(mt::NurbsCurve::new(c))) -} - -fn surface(knots: (mt::KnotVec, mt::KnotVec), rows: Vec>) -> Result { - let s = mt::BSplineSurface::try_new(knots, rows).map_err(fail)?; - Ok(mt::Surface::NurbsSurface(mt::NurbsSurface::new(s))) -} - /// One face from a clamped tensor surface, bounded by its four boundary curves. fn tensor_face(n: &NurbsSurface) -> Result { let [nu, nv] = n.count(); @@ -69,8 +56,7 @@ fn tensor_face(n: &NurbsSurface) -> Result { let top = edge(v3, v2, nurbs_curve(ku.clone(), &row(nv - 1))?)?; let left = edge(v0, v3, nurbs_curve(kv.clone(), &column(0))?)?; let wire: mt::Wire = vec![bottom, right, top.inverse(), left.inverse()].into(); - let rows: Vec> = (0..nu).map(|i| (0..nv).map(|j| h4(n.point(i, j))).collect()).collect(); - mt::Face::try_new(vec![wire], surface((ku, kv), rows)?).map_err(fail) + mt::Face::try_new(vec![wire], crate::spline_conv::surface(n)?).map_err(fail) } /// A parameter point, exact (knot values are shared, so equal points have equal bits). @@ -148,7 +134,7 @@ fn patch_faces(patches: &[BezierPatch]) -> Result> { } let [ps, pt] = p.degree; let rows: Vec> = (0..=ps).map(|i| (0..=pt).map(|j| h4(p.point(i, j))).collect()).collect(); - let surf = surface((mt::KnotVec::bezier_knot(ps), mt::KnotVec::bezier_knot(pt)), rows)?; + let surf = nurbs_surface((mt::KnotVec::bezier_knot(ps), mt::KnotVec::bezier_knot(pt)), rows)?; faces.push(mt::Face::try_new(vec![wire.into()], surf).map_err(fail)?); } Ok(faces) diff --git a/crates/ui-egui/src/preselect_tests.rs b/crates/ui-egui/src/preselect_tests.rs index 241ff237..b301e903 100644 --- a/crates/ui-egui/src/preselect_tests.rs +++ b/crates/ui-egui/src/preselect_tests.rs @@ -709,6 +709,31 @@ fn surface_dialogs_make_their_commands() { assert!((v - 30.0 * 20.0 * 2.0).abs() < 1.0, "{v}"); } +#[test] +fn patch_dialog_fills_a_loop_of_edges_that_is_not_planar() { + // A box with its top cut away by a sideways cylinder: its rim is a closed loop of four + // curved edges, not in a plane. Picking them in the Patch dialog fills it with one face. + let mut app = app_with(json!([ + {"command": "solid.box", "params": {"length": 28, "width": 28, "height": 20, "corner": [-14, -14, 0], "body_name": "Saddle"}}, + {"command": "solid.cylinder", "params": {"radius": 40, "height": 120, "base": [0, -60, 50], "axis": "Y", "operation": "cut", "targets": ["Saddle"]}} + ])); + let z = |x: f64| 50.0 - (1600.0 - x * x).sqrt(); + let rim = [Vec3::new(14.0, 0.0, z(14.0)), Vec3::new(0.0, 14.0, z(0.0)), Vec3::new(-14.0, 0.0, z(-14.0)), Vec3::new(0.0, -14.0, z(0.0))]; + let mut d = start(&mut app, vec![], "surface.patch"); + d.inputs[0].items = rim.iter().enumerate().map(|(i, p)| Sel::Edge { body: "Saddle".into(), index: i, point: *p }).collect(); + let cmds = apply_commands(&app, &d).unwrap(); + assert_eq!(cmds.len(), 1); + assert_eq!(cmds[0].0, "surface.patch"); + app.session.execute(&cmds[0].0, &cmds[0].1).unwrap(); + let st = app.session.model.state(); + let fill = st.bodies.iter().find(|b| b.name != "Saddle").unwrap(); + assert!(fill.body.is_surface() && fill.body.face_count() == 1); + // One smooth face over the rim: its area is close to the cylinder's inside the rim. + let area = fill.mesh().measure().area; + let want = 28.0 * 40.0 * 2.0 * (14.0f64 / 40.0).asin(); + assert!((area - want).abs() < want * 1e-2, "{area} vs {want}"); +} + #[test] fn split_body_by_a_curved_face_sends_it_as_a_tool() { let mut app = app_with(json!([