From a04ef4ca595e2213c6cfa811e214b00b30b38a62 Mon Sep 17 00:00:00 2001 From: Iraklis Kourtis Date: Mon, 12 Oct 2026 13:46:22 -0400 Subject: [PATCH] T-spline core library and sheet bodies (#250, step 1) (#337) * Add B-spline basis, knot insertion and NURBS surfaces to geom Freeform faces need local refinement and gap-free joins between patches with different knot vectors, which tensor B-spline grids cannot give; T-splines can, and the first user is an N-sided Boundary Fill of non-planar loops. This lays the univariate groundwork they rest on, in the L0 geometry crate so it stays independent of the B-rep kernel: - Cox-de Boor basis of one function from its local knot vector, on half-open spans with a left-limit side for the domain end, with derivatives (NURBS Book eq. 2.5, 2.9) - knot insertion into a single basis function (Sederberg et al. 2004 section 3) and Bezier extraction of one function by repeated insertion - Boehm knot insertion for curves (NURBS Book A5.1), de Casteljau evaluation and sub-segments - a rational tensor NURBS surface with evaluation, knot insertion and clamping, and a rational Bezier patch type All input is checked (finite, monotone knots, positive weights, size caps) and returns SplineError; nothing indexes unchecked. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01J53sR6gEaVt56soFDpTck8 * Add T-spline surfaces: T-mesh, evaluation, AS check, refinement, extraction Tensor B-spline faces refine only by whole rows and columns and cannot join patches with different knot vectors without gaps. T-splines allow T-junctions in the control net, which fixes both; the first user is an N-sided Boundary Fill. This adds the maths only (no command, no change in behaviour): - T-mesh in index space with knot values per line, validated (outer rectangle, rectangular faces, no empty lines, size caps); local knot vectors read off rays from each vertex, rays leaving the mesh repeat the boundary (Sederberg et al. 2003 section 3) - rational T-spline evaluation and first derivatives for odd bi-degree, with a complete (p-1)/2-line frame and the domain inside it - T-junction extensions and the analysis-suitability check (Scott et al. 2012 section 2) - local refinement that keeps the surface exactly: requested knot lines are cut in, then every original blending function is carried as a record and split (violation 1) or the mesh gains edges (violations 2 and 3) until each record is a blending function of the refined mesh (Sederberg et al. 2004 section 3); analysis-suitable input stays analysis-suitable by carrying conflicting T-junctions across faces - NURBS <-> T-spline conversion for meshes without T-junctions, and Bezier extraction over the rectangles no blending-function knot line crosses Measured: regular T-spline vs NURBS < 1e-12, round trip exact, partition of unity with T-junctions < 1e-14, refinement deviation <= 1.4e-14 over 40 random refinements, Bezier patches vs surface < 1e-14. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01J53sR6gEaVt56soFDpTck8 * Convert T-splines to sewn NURBS sheet bodies in the kernel A T-spline is only useful to the modeller once it is a B-rep face set. tspline_sheet builds an open shell: a mesh without T-junctions becomes one clamped NURBS face; otherwise each Bezier element is a face, and sides are split at the T-points where neighbours' corners land, so all faces share vertices and edges and the sheet is one oriented shell (V - E + F = 1). Every truck call runs under guard. No command yet; Boundary Fill will be the first caller. Tests: one face for a regular surface, sewn Bezier faces with T-junctions, tessellated points on the surface to ~1e-9 mm, and hostile control nets (coincident, 1e300, needle-thin) never panic. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01J53sR6gEaVt56soFDpTck8 --------- Co-authored-by: Claude Opus 5.5 --- Cargo.lock | 1 + crates/geom/Cargo.toml | 1 + crates/geom/src/bspline.rs | 812 +++++++++++++++++++++++++++++ crates/geom/src/lib.rs | 8 +- crates/geom/src/tspline/extract.rs | 148 ++++++ crates/geom/src/tspline/mesh.rs | 442 ++++++++++++++++ crates/geom/src/tspline/mod.rs | 219 ++++++++ crates/geom/src/tspline/refine.rs | 317 +++++++++++ crates/geom/src/tspline/tests.rs | 350 +++++++++++++ crates/kernel/src/lib.rs | 2 + crates/kernel/src/surfaces.rs | 2 +- crates/kernel/src/tspline.rs | 267 ++++++++++ 12 files changed, 2566 insertions(+), 3 deletions(-) create mode 100644 crates/geom/src/bspline.rs create mode 100644 crates/geom/src/tspline/extract.rs create mode 100644 crates/geom/src/tspline/mesh.rs create mode 100644 crates/geom/src/tspline/mod.rs create mode 100644 crates/geom/src/tspline/refine.rs create mode 100644 crates/geom/src/tspline/tests.rs create mode 100644 crates/kernel/src/tspline.rs diff --git a/Cargo.lock b/Cargo.lock index 187fdfa7..898a3fba 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3351,6 +3351,7 @@ version = "0.4.0" dependencies = [ "serde", "serde_json", + "thiserror 2.0.21", ] [[package]] diff --git a/crates/geom/Cargo.toml b/crates/geom/Cargo.toml index 2f5c6b18..0c1976ba 100644 --- a/crates/geom/Cargo.toml +++ b/crates/geom/Cargo.toml @@ -10,6 +10,7 @@ workspace = true [dependencies] serde = { workspace = true } +thiserror = { workspace = true } [dev-dependencies] serde_json = { workspace = true } diff --git a/crates/geom/src/bspline.rs b/crates/geom/src/bspline.rs new file mode 100644 index 00000000..a8e0f968 --- /dev/null +++ b/crates/geom/src/bspline.rs @@ -0,0 +1,812 @@ +//! B-spline machinery shared by NURBS and T-spline surfaces: homogeneous points, the +//! univariate basis (Cox–de Boor on half-open intervals, with derivatives), Boehm knot +//! insertion, knot insertion into a single basis function, Bézier extraction, and a tensor +//! NURBS surface. +//! +//! References: Piegl & Tiller, *The NURBS Book* (2nd ed.) — §2.2 (basis definition, eq. 2.5), +//! §2.3 (derivatives, eq. 2.9), A2.1 (span search), §5.2 (knot insertion, eq. 5.10 / A5.1), +//! §1.4 (de Casteljau); Sederberg, Cardon, Finnigan, North, Zheng, Lyche 2004, "T-spline +//! simplification and local refinement", §3 (knot insertion into one blending function). +//! +//! Input is hostile: every constructor checks finiteness and sizes and returns +//! [`SplineError`] instead of panicking; no function here indexes with an unchecked number. + +use std::ops::{Add, AddAssign, Mul, Sub}; + +use crate::Vec3; + +/// Highest polynomial degree accepted in either direction. +pub const MAX_SPLINE_DEGREE: usize = 7; +/// Most knots accepted in one knot vector. +pub const MAX_KNOTS: usize = 1 << 14; +/// Most control points accepted in one tensor surface. +pub const MAX_SURFACE_POINTS: usize = 1 << 20; + +/// Why a spline could not be built or used. +#[derive(Debug, Clone, PartialEq, thiserror::Error)] +pub enum SplineError { + #[error("invalid spline: {0}")] + Invalid(String), + #[error("not supported yet: {0}")] + Unsupported(String), + #[error("spline too large: {0}")] + TooLarge(String), +} + +pub type SplineResult = Result; + +pub(crate) fn invalid(msg: impl Into) -> SplineResult { + Err(SplineError::Invalid(msg.into())) +} + +/// A homogeneous control point `(x·w, y·w, z·w, w)`. +#[derive(Clone, Copy, Debug, Default, PartialEq)] +pub struct HPoint { + pub x: f64, + pub y: f64, + pub z: f64, + pub w: f64, +} + +impl HPoint { + pub const ZERO: HPoint = HPoint { x: 0.0, y: 0.0, z: 0.0, w: 0.0 }; + + /// The point `p` with weight `w`. + pub fn weighted(p: Vec3, w: f64) -> HPoint { + HPoint { x: p.x * w, y: p.y * w, z: p.z * w, w } + } + /// The Cartesian point, if the weight is usable. + pub fn project(self) -> Option { + let p = Vec3::new(self.x / self.w, self.y / self.w, self.z / self.w); + (self.w.abs() > f64::MIN_POSITIVE && p.is_finite()).then_some(p) + } + /// The weighted coordinates `(x·w, y·w, z·w)`. + pub fn xyz(self) -> Vec3 { + Vec3::new(self.x, self.y, self.z) + } + pub fn is_finite(self) -> bool { + self.x.is_finite() && self.y.is_finite() && self.z.is_finite() && self.w.is_finite() + } +} + +impl Add for HPoint { + type Output = HPoint; + fn add(self, o: HPoint) -> HPoint { + HPoint { x: self.x + o.x, y: self.y + o.y, z: self.z + o.z, w: self.w + o.w } + } +} +impl Sub for HPoint { + type Output = HPoint; + fn sub(self, o: HPoint) -> HPoint { + HPoint { x: self.x - o.x, y: self.y - o.y, z: self.z - o.z, w: self.w - o.w } + } +} +impl Mul for HPoint { + type Output = HPoint; + fn mul(self, k: f64) -> HPoint { + HPoint { x: self.x * k, y: self.y * k, z: self.z * k, w: self.w * k } + } +} +impl AddAssign for HPoint { + fn add_assign(&mut self, o: HPoint) { + *self = *self + o; + } +} + +/// A parameter direction: along s (u, the first) or along t (v, the second). +#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)] +pub enum Dir { + S, + T, +} + +impl Dir { + /// `s` for [`Dir::S`], `t` for [`Dir::T`]. + pub fn pick(self, s: X, t: X) -> X { + match self { + Dir::S => s, + Dir::T => t, + } + } +} + +/// Homogeneous control points from finite points and positive weights, one each. +pub fn weighted_points(points: &[Vec3], weights: &[f64]) -> SplineResult> { + if points.len() != weights.len() { + return invalid("one weight per control point"); + } + if points.iter().any(|p| !p.is_finite()) { + return invalid("a control point is not finite"); + } + if weights.iter().any(|w| !(w.is_finite() && *w > 0.0)) { + return invalid("weights must be positive numbers"); + } + Ok(points.iter().zip(weights).map(|(p, w)| HPoint::weighted(*p, *w)).collect()) +} + +/// Which one-sided limit a basis evaluation takes at a knot. Basis functions live on half-open +/// spans `[uᵢ, uᵢ₊₁)` ([`Side::Right`]); at the right end of a domain the left limit on +/// `(uᵢ, uᵢ₊₁]` ([`Side::Left`]) keeps the closing span. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Side { + Right, + Left, +} + +/// Check a knot vector: finite, non-decreasing, not too long. +pub fn check_knots(knots: &[f64]) -> SplineResult<()> { + if knots.len() > MAX_KNOTS { + return Err(SplineError::TooLarge(format!("{} knots (at most {MAX_KNOTS})", knots.len()))); + } + if knots.iter().any(|k| !k.is_finite()) { + return invalid("a knot is not a finite number"); + } + if knots.windows(2).any(|w| matches!(w, [a, b] if b < a)) { + return invalid("knots must not decrease"); + } + Ok(()) +} + +/// One B-spline basis function of degree `p = knots.len() − 2` given by its local knot vector +/// `[ξ₀ … ξₚ₊₁]`, and its derivatives up to order `n`: element `k` is the `k`-th derivative at +/// `u`. Cox–de Boor recursion (NURBS Book eq. 2.5) with the 0/0 = 0 convention, derivatives by +/// eq. 2.9 applied repeatedly (`dᵏNᵢ,ₚ = p·(dᵏ⁻¹Nᵢ,ₚ₋₁/(uᵢ₊ₚ−uᵢ) − dᵏ⁻¹Nᵢ₊₁,ₚ₋₁/(uᵢ₊ₚ₊₁−uᵢ₊₁))`). +/// 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 Some(p) = knots.len().checked_sub(2) else { return out }; + if p > MAX_SPLINE_DEGREE || !u.is_finite() { + return out; + } + let kn = |i: usize| knots.get(i).copied().unwrap_or(0.0); + let ratio = |num: f64, den: f64| if den > 0.0 { num / den } else { 0.0 }; + // 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) { + let (a, b) = (kn(j), kn(j + 1)); + let inside = match side { + Side::Right => a <= u && u < b, + Side::Left => a < u && u <= b, + }; + *slot = if inside { 1.0 } else { 0.0 }; + } + for q in 1..=p { + for k in (0..=levels).rev() { + for j in 0..=p - q { + let (d1, d2) = (kn(j + q) - kn(j), kn(j + q + 1) - kn(j + 1)); + let at = |k: usize, j: usize| table.get(k).and_then(|r| r.get(j)).copied().unwrap_or(0.0); + let v = if k == 0 { + ratio(u - kn(j), d1) * at(0, j) + ratio(kn(j + q + 1) - u, d2) * at(0, j + 1) + } else { + q as f64 * (ratio(at(k - 1, j), d1) - ratio(at(k - 1, j + 1), d2)) + }; + if let Some(slot) = table.get_mut(k).and_then(|r| r.get_mut(j)) { + *slot = v; + } + } + } + } + for (k, o) in out.iter_mut().enumerate().take(levels + 1) { + *o = table.get(k).map_or(0.0, |r| r[0]); + } + out +} + +/// The coefficients `(c₁, c₂)` of knot insertion into one basis function: inserting `k` +/// (with `ξ₀ ≤ k ≤ ξₚ₊₁`) into `N[ξ₀…ξₚ₊₁]` gives `c₁·N[first p+2 merged knots] + +/// c₂·N[last p+2 merged knots]` (Sederberg et al. 2004, §3; the single-function case of the +/// NURBS Book eq. 5.10). `None` when `k` is outside the support or the knots are too few. +pub fn split_basis(knots: &[f64], k: f64) -> Option<(f64, f64)> { + let p = knots.len().checked_sub(2)?; + let (x0, x1, xp, xl) = (*knots.first()?, *knots.get(1)?, *knots.get(p)?, *knots.get(p + 1)?); + if !(k >= x0 && k <= xl) { + return None; + } + let c1 = if k >= xp { 1.0 } else { (k - x0) / (xp - x0) }; + let c2 = if k <= x1 { 1.0 } else { (xl - k) / (xl - x1) }; + Some((c1, c2)) +} + +/// The two knot vectors of a split basis function (see [`split_basis`]): `k` merged into +/// `knots` in order, then the first and the last `knots.len()` entries. +pub fn split_knots(knots: &[T], k: T) -> [Vec; 2] { + let at = knots.iter().position(|&x| x > k).unwrap_or(knots.len()); + let mut merged = knots.to_vec(); + merged.insert(at, k); + let n = knots.len(); + [merged.get(..n).map(<[T]>::to_vec).unwrap_or_default(), merged.get(1..).map(<[T]>::to_vec).unwrap_or_default()] +} + +/// The Bernstein coefficients of one basis function `N[knots]` on `[a, b]`: `N(u) = +/// Σ cᵢ Bᵢ,ₚ((u−a)/(b−a))` there. `[a, b]` must not contain a knot in its interior (the +/// function is one polynomial on it). Found by inserting `a` and `b` into the function until +/// every surviving piece is a Bernstein polynomial `N[a…a, b…b]` (Sederberg et al. 2004, §3; +/// NURBS Book §5.4, decomposition into Bézier segments). +pub fn basis_bezier(knots: &[f64], a: f64, b: f64) -> SplineResult> { + let Some(p) = knots.len().checked_sub(2) else { return invalid("a basis function needs at least two knots") }; + if p > MAX_SPLINE_DEGREE { + return Err(SplineError::Unsupported(format!("degree {p} (at most {MAX_SPLINE_DEGREE})"))); + } + check_knots(knots)?; + if !(a.is_finite() && b.is_finite() && a < b) { + return invalid("a Bézier interval must be finite and non-empty"); + } + if knots.iter().any(|&x| x > a && x < b) { + return invalid("a knot lies inside the Bézier interval"); + } + let mut out = vec![0.0; p + 1]; + let mut work = vec![(knots.to_vec(), 1.0)]; + let mut steps = 0usize; + while let Some((kn, c)) = work.pop() { + steps += 1; + if steps > 1 << 12 { + return Err(SplineError::TooLarge("Bézier extraction did not converge".into())); + } + let (Some(&lo), Some(&hi)) = (kn.first(), kn.last()) else { continue }; + if hi <= a || lo >= b || c == 0.0 { + continue; + } + let insert = if kn.iter().any(|&x| x < a) { + Some(a) + } else if kn.iter().any(|&x| x > b) { + Some(b) + } else { + None + }; + match insert { + Some(k) => { + let (c1, c2) = split_basis(&kn, k).ok_or_else(|| SplineError::Invalid("knot insertion outside the support".into()))?; + let [first, second] = split_knots(&kn, k); + work.push((first, c * c1)); + work.push((second, c * c2)); + } + None => { + // Knots a×(p+1−i), b×(i+1): the Bernstein polynomial Bᵢ,ₚ. + let nb = kn.iter().filter(|&&x| x == b).count(); + if let Some(slot) = nb.checked_sub(1).and_then(|i| out.get_mut(i)) { + *slot += c; + } + } + } + } + Ok(out) +} + +/// The span index `k` with `knots[k] ≤ u < knots[k+1]` among the valid spans `p ≤ k < n` +/// (`n` control points), the last span closed at the domain end (NURBS Book A2.1). +pub fn find_span(knots: &[f64], p: usize, n: usize, u: f64) -> Option { + if n <= p || knots.len() != n + p + 1 || !u.is_finite() { + return None; + } + let (lo, hi) = (*knots.get(p)?, *knots.get(n)?); + 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; + } + } + Some(k) +} + +/// Insert the knot `u` once into a B-spline with degree `p`, knots and homogeneous control +/// points (Boehm; NURBS Book eq. 5.10 / A5.1). `u` must lie in the domain `[knots[p], +/// knots[n]]`; the curve is unchanged. +pub fn insert_knot(p: usize, knots: &[f64], ctrl: &[HPoint], u: f64) -> SplineResult<(Vec, Vec)> { + let n = ctrl.len(); + if n <= p || knots.len() != n + p + 1 { + return invalid("knot and control point counts do not match the degree"); + } + if knots.len() >= MAX_KNOTS { + return Err(SplineError::TooLarge("too many knots".into())); + } + check_knots(knots)?; + let k = find_span(knots, p, n, u).ok_or_else(|| SplineError::Invalid(format!("knot {u} is outside the domain")))?; + let kn = |i: usize| knots.get(i).copied().unwrap_or(0.0); + let pt = |i: usize| ctrl.get(i).copied().unwrap_or(HPoint::ZERO); + let mut q = Vec::with_capacity(n + 1); + for i in 0..=n { + q.push(if i + p <= k { + pt(i) + } else if i > k { + pt(i - 1) + } else { + let den = kn(i + p) - kn(i); + let alpha = if den > 0.0 { (u - kn(i)) / den } else { 1.0 }; + pt(i) * alpha + pt(i - 1) * (1.0 - alpha) + }); + } + let mut nk = knots.to_vec(); + nk.insert(k + 1, u); + Ok((nk, q)) +} + +/// 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; + } + c * x.powi(i as i32) * (1.0 - x).powi((p - i) as i32) +} + +/// A Bézier curve's point at `x ∈ [0, 1]` (de Casteljau, NURBS Book §1.4). +pub fn bezier_point(ctrl: &[HPoint], x: f64) -> HPoint { + subdivide(ctrl, x).0.last().copied().unwrap_or(HPoint::ZERO) +} + +/// The control points of the piece of a Bézier curve over `[a, b] ⊂ [0, 1]` (two de Casteljau +/// subdivisions). +pub fn bezier_segment(ctrl: &[HPoint], a: f64, b: f64) -> SplineResult> { + if !(a.is_finite() && b.is_finite() && 0.0 <= a && a < b && b <= 1.0) { + return invalid("a Bézier segment must satisfy 0 ≤ a < b ≤ 1"); + } + // Left part [0, b], then its right part from a/b. + let left = subdivide(ctrl, b).0; + Ok(subdivide(&left, a / b).1) +} + +/// Split a Bézier curve at `x` (de Casteljau): control points of `[0, x]` and `[x, 1]`. +fn subdivide(ctrl: &[HPoint], x: f64) -> (Vec, Vec) { + let mut left = Vec::with_capacity(ctrl.len()); + let mut right = Vec::with_capacity(ctrl.len()); + let mut pts = ctrl.to_vec(); + while let (Some(&f), Some(&l)) = (pts.first(), pts.last()) { + left.push(f); + right.push(l); + pts = pts.windows(2).map(|w| if let [a, b] = w { *a * (1.0 - x) + *b * x } else { HPoint::ZERO }).collect(); + } + right.reverse(); + (left, right) +} + +/// A rational tensor Bézier patch over the parameter rectangle `s × t`: control point `(i, j)` +/// (i along s, j along t) at `ctrl[j·(degree[0]+1) + i]`. +#[derive(Clone, Debug, PartialEq)] +pub struct BezierPatch { + pub s: [f64; 2], + pub t: [f64; 2], + pub degree: [usize; 2], + pub ctrl: Vec, +} + +impl BezierPatch { + /// Control point `(i, j)`. + pub fn point(&self, i: usize, j: usize) -> HPoint { + self.ctrl.get(j * (self.degree[0] + 1) + i).copied().unwrap_or(HPoint::ZERO) + } + /// Homogeneous point at global parameters `(s, t)`. + pub fn eval_h(&self, s: f64, t: f64) -> HPoint { + let x = (s - self.s[0]) / (self.s[1] - self.s[0]); + let y = (t - self.t[0]) / (self.t[1] - self.t[0]); + let mut acc = HPoint::ZERO; + for j in 0..=self.degree[1] { + let bj = bernstein(self.degree[1], j, y); + for i in 0..=self.degree[0] { + acc += self.point(i, j) * (bernstein(self.degree[0], i, x) * bj); + } + } + acc + } + /// Point at global parameters `(s, t)`. + pub fn eval(&self, s: f64, t: f64) -> Option { + self.eval_h(s, t).project() + } + /// The control points of one side as a Bézier curve: `0` t = t₀ (s rising), `1` s = s₁ + /// (t rising), `2` t = t₁ (s rising), `3` s = s₀ (t rising). + pub fn side(&self, side: usize) -> Vec { + let [ps, pt] = self.degree; + match side { + 0 => (0..=ps).map(|i| self.point(i, 0)).collect(), + 1 => (0..=pt).map(|j| self.point(ps, j)).collect(), + 2 => (0..=ps).map(|i| self.point(i, pt)).collect(), + _ => (0..=pt).map(|j| self.point(0, j)).collect(), + } + } +} + +/// A rational tensor-product B-spline (NURBS) surface. Control point `(i, j)` (i along u, j +/// along v) is `ctrl[j·count_u + i]`; the domain is `[knots_u[p], knots_u[n]] × …`. +#[derive(Clone, Debug, PartialEq)] +pub struct NurbsSurface { + degree: [usize; 2], + knots: [Vec; 2], + count: [usize; 2], + ctrl: Vec, +} + +impl NurbsSurface { + /// A surface from Cartesian points and positive weights (one per point, u fastest). + pub fn new(degree: [usize; 2], knots: [Vec; 2], points: &[Vec3], weights: &[f64]) -> SplineResult { + Self::from_homogeneous(degree, knots, weighted_points(points, weights)?) + } + + /// A surface from homogeneous control points (u fastest). + pub fn from_homogeneous(degree: [usize; 2], knots: [Vec; 2], ctrl: Vec) -> SplineResult { + let mut count = [0; 2]; + for d in 0..2 { + let (p, kn) = (degree[d], &knots[d]); + if p == 0 || p > MAX_SPLINE_DEGREE { + return Err(SplineError::Unsupported(format!("degree {p} (1 to {MAX_SPLINE_DEGREE})"))); + } + check_knots(kn)?; + let n = kn.len().checked_sub(p + 1).filter(|&n| n > p).ok_or_else(|| SplineError::Invalid("too few knots for the degree".into()))?; + if !(kn.get(p) < kn.get(n)) { + return invalid("the surface's domain is empty"); + } + count[d] = n; + } + if count[0].checked_mul(count[1]).is_none_or(|c| c > MAX_SURFACE_POINTS) { + return Err(SplineError::TooLarge(format!("more than {MAX_SURFACE_POINTS} control points"))); + } + if ctrl.len() != count[0] * count[1] { + return invalid(format!("{} control points, expected {}", ctrl.len(), count[0] * count[1])); + } + if ctrl.iter().any(|h| !h.is_finite() || !(h.w > 0.0)) { + return invalid("control points must be finite with positive weights"); + } + Ok(NurbsSurface { degree, knots, count, ctrl }) + } + + pub fn degree(&self) -> [usize; 2] { + self.degree + } + pub fn knots(&self) -> &[Vec; 2] { + &self.knots + } + /// Control point counts along u and v. + pub fn count(&self) -> [usize; 2] { + self.count + } + pub fn homogeneous(&self) -> &[HPoint] { + &self.ctrl + } + /// Control point `(i, j)`. + pub fn point(&self, i: usize, j: usize) -> HPoint { + self.ctrl.get(j * self.count[0] + i).copied().unwrap_or(HPoint::ZERO) + } + + /// The parameter domain `[[u₀, u₁], [v₀, v₁]]`. + pub fn domain(&self) -> [[f64; 2]; 2] { + let d = |k: usize| { + let (p, n, kn) = (self.degree[k], self.count[k], &self.knots[k]); + [kn.get(p).copied().unwrap_or(0.0), kn.get(n).copied().unwrap_or(0.0)] + }; + [d(0), d(1)] + } + + /// The point at `(u, v)`, which must lie in the domain. + pub fn eval(&self, u: f64, v: f64) -> SplineResult { + let dom = self.domain(); + let mut spans = [0; 2]; + let mut basis: [Vec; 2] = [Vec::new(), Vec::new()]; + for (d, x) in [u, v].into_iter().enumerate() { + let (p, n, kn) = (self.degree[d], self.count[d], &self.knots[d]); + let span = find_span(kn, p, n, x).ok_or_else(|| SplineError::Invalid(format!("({u}, {v}) is outside the domain")))?; + let side = if x >= dom[d][1] { Side::Left } else { Side::Right }; + spans[d] = span; + basis[d] = (span - p..=span).map(|i| kn.get(i..i + p + 2).map_or(0.0, |k| local_basis(k, x, side, 0)[0])).collect(); + } + let mut acc = HPoint::ZERO; + for (b, nj) in basis[1].iter().enumerate() { + for (a, ni) in basis[0].iter().enumerate() { + acc += self.point(spans[0] - self.degree[0] + a, spans[1] - self.degree[1] + b) * (ni * nj); + } + } + acc.project().ok_or_else(|| SplineError::Invalid("the surface's weights vanish there".into())) + } + + /// The same surface with the knot `x` inserted once in direction `dir` ([`Dir::S`] = u). + pub fn insert_knot(&self, dir: Dir, x: f64) -> SplineResult { + let d = dir.pick(0, 1); + let (p, kn) = (self.degree[d], &self.knots[d]); + let [nu, nv] = self.count; + let (lines, len) = if d == 0 { (nv, nu) } else { (nu, nv) }; + let at = |line: usize, k: usize| if d == 0 { self.point(k, line) } else { self.point(line, k) }; + let mut new_knots = Vec::new(); + let mut rows: Vec> = Vec::with_capacity(lines); + for line in 0..lines { + let row: Vec = (0..len).map(|k| at(line, k)).collect(); + let (nk, nr) = insert_knot(p, kn, &row, x)?; + new_knots = nk; + rows.push(nr); + } + let mut knots = self.knots.clone(); + knots[d] = new_knots; + let new_len = len + 1; + let mut ctrl = Vec::with_capacity(lines * new_len); + if d == 0 { + for row in &rows { + ctrl.extend_from_slice(row); + } + } else { + for j in 0..new_len { + for row in &rows { + ctrl.push(row.get(j).copied().unwrap_or(HPoint::ZERO)); + } + } + } + NurbsSurface::from_homogeneous(self.degree, knots, ctrl) + } + + /// 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). + pub fn clamped(&self) -> SplineResult { + let mut s = self.clone(); + for dir in [Dir::S, Dir::T] { + let d = dir.pick(0, 1); + 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)?; + } + } + 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)] +pub(crate) mod tests { + use super::*; + + /// A small deterministic generator for test data. + pub(crate) struct Lcg(pub(crate) u64); + impl Lcg { + pub(crate) fn unit(&mut self) -> f64 { + self.0 = self.0.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + (self.0 >> 11) as f64 / (1u64 << 53) as f64 + } + } + + const KNOTS: [f64; 12] = [0.0, 0.0, 0.0, 0.0, 0.3, 0.3, 0.5, 0.9, 1.0, 1.0, 1.0, 1.0]; + + fn samples() -> Vec { + (0..=200).map(|k| k as f64 / 200.0).collect() + } + + #[test] + fn basis_partition_of_unity_non_negative_compact_support() { + let p = 3; + let n = KNOTS.len() - p - 1; + for u in samples() { + let side = if u >= 1.0 { Side::Left } else { Side::Right }; + let mut sum = 0.0; + for i in 0..n { + let local = &KNOTS[i..i + p + 2]; + let b = local_basis(local, u, side, 0)[0]; + assert!(b >= -1e-15, "N{i}({u}) = {b}"); + if u < local[0] || u > local[p + 1] { + assert_eq!(b, 0.0, "N{i} is not zero outside its support at {u}"); + } + sum += b; + } + assert!((sum - 1.0).abs() < 1e-14, "Σ N({u}) = {sum}"); + } + } + + #[test] + fn basis_derivatives_match_finite_differences() { + let p = 3; + let h = 1e-6; + let mut worst = 0.0f64; + for i in 0..KNOTS.len() - p - 1 { + let local = &KNOTS[i..i + p + 2]; + for u in samples() { + if KNOTS.iter().any(|k| (k - u).abs() < 2.0 * h) { + continue; + } + let d = local_basis(local, u, Side::Right, 2); + let f = |x: f64| local_basis(local, x, Side::Right, 0)[0]; + let fd1 = (f(u + h) - f(u - h)) / (2.0 * h); + let df = |x: f64| local_basis(local, x, Side::Right, 1)[1]; + let fd2 = (df(u + h) - df(u - h)) / (2.0 * h); + worst = worst.max((d[1] - fd1).abs()).max((d[2] - fd2).abs() / 100.0); + } + } + assert!(worst < 1e-5, "derivative error {worst}"); + } + + #[test] + fn single_function_knot_insertion_is_exact() { + let mut rng = Lcg(7); + for _ in 0..200 { + let mut kn: Vec = (0..5).map(|_| (rng.unit() * 8.0).round() / 8.0).collect(); + kn.sort_by(f64::total_cmp); + if kn[0] == kn[4] { + continue; + } + let k = kn[0] + (kn[4] - kn[0]) * (rng.unit() * 8.0).round() / 8.0; + let (c1, c2) = split_basis(&kn, k).unwrap(); + let mut merged = kn.clone(); + let at = merged.iter().position(|&x| x > k).unwrap_or(5); + merged.insert(at, k); + for u in samples() { + let u = kn[0] - 0.1 + u * (kn[4] - kn[0] + 0.2); + let f = |k: &[f64]| local_basis(k, u, Side::Right, 0)[0]; + let lhs = f(&kn); + let rhs = c1 * f(&merged[..5]) + c2 * f(&merged[1..]); + assert!((lhs - rhs).abs() < 1e-14, "{kn:?} insert {k} at {u}: {lhs} vs {rhs}"); + } + } + } + + #[test] + fn basis_bezier_coefficients_reproduce_the_function() { + let p = 3; + let spans: Vec<(f64, f64)> = KNOTS.windows(2).filter(|w| w[0] < w[1]).map(|w| (w[0], w[1])).collect(); + for i in 0..KNOTS.len() - p - 1 { + let local = &KNOTS[i..i + p + 2]; + for &(a, b) in &spans { + let c = basis_bezier(local, a, b).unwrap(); + for k in 0..=20 { + let x = k as f64 / 20.0; + let u = a + (b - a) * x; + let side = if x >= 1.0 { Side::Left } else { Side::Right }; + let want = local_basis(local, u, side, 0)[0]; + let got: f64 = c.iter().enumerate().map(|(j, cj)| cj * bernstein(p, j, x)).sum(); + assert!((want - got).abs() < 1e-14, "N{i} on [{a}, {b}] at {u}: {want} vs {got}"); + } + } + } + assert!(basis_bezier(&KNOTS[..5], 0.0, 0.5).is_err(), "a knot inside the interval is refused"); + assert!(basis_bezier(&KNOTS[..5], f64::NAN, 0.5).is_err()); + } + + fn curve(p: usize, knots: &[f64], ctrl: &[HPoint], u: f64) -> Vec3 { + let n = ctrl.len(); + let side = if u >= knots[n] { Side::Left } else { Side::Right }; + let mut acc = HPoint::ZERO; + for (i, c) in ctrl.iter().enumerate() { + acc += *c * local_basis(&knots[i..i + p + 2], u, side, 0)[0]; + } + acc.project().unwrap() + } + + #[test] + fn boehm_insertion_keeps_the_curve() { + let p = 3; + let mut rng = Lcg(3); + let ctrl: Vec = (0..8).map(|_| HPoint::weighted(Vec3::new(rng.unit(), rng.unit(), rng.unit()) * 10.0, 0.5 + rng.unit())).collect(); + let (mut kn, mut pts) = (KNOTS.to_vec(), ctrl.clone()); + for u in [0.0, 0.1, 0.3, 0.3, 0.77, 1.0, 1.0] { + (kn, pts) = insert_knot(p, &kn, &pts, u).unwrap(); + } + let mut worst = 0.0f64; + for u in samples() { + worst = worst.max(curve(p, &KNOTS, &ctrl, u).dist(curve(p, &kn, &pts, u))); + } + assert!(worst < 1e-13, "deviation {worst}"); + assert!(insert_knot(p, &KNOTS, &ctrl, 1.5).is_err()); + assert!(insert_knot(p, &KNOTS, &ctrl, f64::NAN).is_err()); + } + + pub(crate) fn random_surface(seed: u64, ku: Vec, kv: Vec, degree: [usize; 2]) -> NurbsSurface { + let mut rng = Lcg(seed); + let n = (ku.len() - degree[0] - 1) * (kv.len() - degree[1] - 1); + let pts: Vec = (0..n).map(|_| Vec3::new(rng.unit(), rng.unit(), rng.unit()) * 10.0).collect(); + let w: Vec = (0..n).map(|_| 0.5 + rng.unit()).collect(); + NurbsSurface::new(degree, [ku, kv], &pts, &w).unwrap() + } + + fn max_deviation(a: &NurbsSurface, b: &NurbsSurface) -> f64 { + let [du, dv] = a.domain(); + let mut worst = 0.0f64; + for i in 0..=30 { + for j in 0..=30 { + let u = du[0] + (du[1] - du[0]) * i as f64 / 30.0; + let v = dv[0] + (dv[1] - dv[0]) * j as f64 / 30.0; + worst = worst.max(a.eval(u, v).unwrap().dist(b.eval(u, v).unwrap())); + } + } + worst + } + + #[test] + fn surface_knot_insertion_and_clamping_keep_the_surface() { + let uniform: Vec = (0..10).map(|k| k as f64).collect(); + let s = random_surface(11, uniform.clone(), (0..9).map(|k| k as f64 * 0.5).collect(), [3, 3]); + let r = s.insert_knot(Dir::S, 4.5).unwrap().insert_knot(Dir::T, 2.25).unwrap(); + assert!(max_deviation(&s, &r) < 1e-12); + let c = s.clamped().unwrap(); + assert_eq!(c.domain(), s.domain()); + assert_eq!(&c.knots()[0][..4], &[3.0; 4]); + assert_eq!(&c.knots()[1][c.knots()[1].len() - 4..], &[2.5; 4]); + let dev = max_deviation(&s, &c); + assert!(dev < 1e-12, "clamping moved the surface by {dev}"); + // A clamped surface interpolates its corner control points. + let corner = c.eval(3.0, 1.5).unwrap(); + assert!(corner.dist(c.point(0, 0).project().unwrap()) < 1e-12); + } + + #[test] + fn bezier_segments_lie_on_the_curve() { + let mut rng = Lcg(5); + let ctrl: Vec = (0..4).map(|_| HPoint::weighted(Vec3::new(rng.unit(), rng.unit(), 0.0), 0.5 + rng.unit())).collect(); + let seg = bezier_segment(&ctrl, 0.2, 0.7).unwrap(); + for k in 0..=10 { + let x = k as f64 / 10.0; + let want = bezier_point(&ctrl, 0.2 + 0.5 * x).project().unwrap(); + assert!(bezier_point(&seg, x).project().unwrap().dist(want) < 1e-14); + } + assert!(bezier_segment(&ctrl, 0.7, 0.2).is_err()); + } + + #[test] + fn hostile_spline_input_is_refused() { + let pts = [Vec3::ZERO; 16]; + let w = [1.0; 16]; + let k = vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]; + assert!(NurbsSurface::new([3, 3], [k.clone(), k.clone()], &pts, &w).is_ok()); + assert!(NurbsSurface::new([3, 3], [vec![0.0, f64::NAN, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0], k.clone()], &pts, &w).is_err()); + assert!(NurbsSurface::new([3, 3], [k.clone(), k.clone()], &pts, &[0.0; 16]).is_err()); + assert!(NurbsSurface::new([3, 3], [k.clone(), k.clone()], &pts[..15], &w[..15]).is_err()); + assert!(NurbsSurface::new([0, 3], [k.clone(), k.clone()], &pts, &w).is_err()); + assert!(NurbsSurface::new([9, 3], [k.clone(), k.clone()], &pts, &w).is_err()); + assert!(NurbsSurface::new([3, 3], [vec![1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0], k.clone()], &pts, &w).is_err()); + assert!(check_knots(&vec![0.0; MAX_KNOTS + 1]).is_err()); + let s = NurbsSurface::new([3, 3], [k.clone(), k], &pts, &w).unwrap(); + assert!(s.eval(2.0, 0.5).is_err()); + assert!(s.eval(f64::NAN, 0.5).is_err()); + assert_eq!(local_basis(&[], 0.5, Side::Right, 1), vec![0.0, 0.0]); + } +} diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index 77991ef6..de9b3b8b 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -1,14 +1,18 @@ -//! Geometry foundation for SolveCraft: 2D/3D vectors, sketch planes, axis-aligned boxes and -//! triangle meshes with their measures (volume, area, centroid). Units are millimetres. +//! Geometry foundation for SolveCraft: 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. //! //! This crate has no knowledge of sketches, documents or the B-rep kernel; every other crate //! builds on it. #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] #![forbid(unsafe_code)] +pub mod bspline; mod mesh; mod plane; mod profile; +pub mod tspline; mod vec; pub use mesh::{Aabb3, InsideIndex, Mesh, MeshMeasure}; diff --git a/crates/geom/src/tspline/extract.rs b/crates/geom/src/tspline/extract.rs new file mode 100644 index 00000000..9a11aa74 --- /dev/null +++ b/crates/geom/src/tspline/extract.rs @@ -0,0 +1,148 @@ +//! Conversions: a tensor NURBS surface as a T-spline and back (a T-mesh without T-junctions is +//! a tensor B-spline; 2003 §3 and §5), and Bézier extraction — the surface as rational Bézier +//! patches over the rectangles where every blending function is one polynomial. + +use super::mesh::{Grid, TMesh}; +use super::{Blend, TSpline, frame}; +use crate::bspline::{BezierPatch, HPoint, NurbsSurface, SplineError, SplineResult, basis_bezier, invalid}; + +/// Most unit-edge marks and coefficient products in one Bézier extraction. +const MAX_EXTRACT_WORK: usize = 1 << 26; + +impl TSpline { + /// A tensor NURBS surface of odd bi-degree as a T-spline on the full grid: column `i` of the + /// mesh is knot `i + (p+1)/2` of the knot vector, so the blending function at a vertex is + /// the B-spline basis function of that control point (rays that leave the mesh repeat the + /// boundary line, 2003 §3). Knot vectors whose first and last `(p+1)/2 + 1` knots are not + /// equal are clamped first (same surface on the same domain). + pub fn from_nurbs(n: &NurbsSurface) -> SplineResult { + let degree = n.degree(); + if degree.iter().any(|p| p % 2 == 0) { + return Err(SplineError::Unsupported("a NURBS surface of even degree as a T-spline".into())); + } + let ends_ok = |k: &[f64], p: usize| { + let h = p.div_ceil(2); + let same = |w: Option<&[f64]>| w.is_some_and(|w| w.windows(2).all(|x| matches!(x, [a, b] if a == b))); + same(k.get(..=h)) && same(k.len().checked_sub(h + 1).and_then(|s| k.get(s..))) + }; + let n = if (0..2).all(|d| ends_ok(&n.knots()[d], degree[d])) { n.clone() } else { n.clamped()? }; + let lines = |d: usize| -> Vec { + let h = degree[d].div_ceil(2); + let k = &n.knots()[d]; + k.get(h..k.len().saturating_sub(h)).map(<[f64]>::to_vec).unwrap_or_default() + }; + let mesh = TMesh::tensor(lines(0), lines(1))?; + TSpline::from_homogeneous(mesh, degree, n.homogeneous().to_vec()) + } + + /// The T-spline as a tensor NURBS surface, when its mesh has no T-junctions. + pub fn to_nurbs(&self) -> SplineResult { + if !self.mesh.is_regular() { + return Err(SplineError::Unsupported("a T-mesh with T-junctions as one tensor surface".into())); + } + let pad = |k: &[f64], p: usize| -> Vec { + let h = p.div_ceil(2); + let (first, last) = (k.first().copied().unwrap_or(0.0), k.last().copied().unwrap_or(0.0)); + std::iter::repeat_n(first, h).chain(k.iter().copied()).chain(std::iter::repeat_n(last, h)).collect() + }; + let g = &self.mesh.grid; + NurbsSurface::from_homogeneous(self.degree, [pad(&g.s, self.degree[0]), pad(&g.t, self.degree[1])], self.ctrl.clone()) + } + + /// The surface as rational Bézier patches over rectangles of the domain, inside each of + /// which every blending function is one polynomial: the knot lines of all blending + /// functions (each function's s-knots across its t-support and t-knots across its + /// s-support) are marked on the index grid and the domain's cells are gathered into + /// rectangles that no marked line crosses. Each function's Bernstein coefficients on a + /// rectangle come from knot insertion ([`basis_bezier`]). Patches tile the domain and meet + /// at T-points where their corners land on a neighbour's side. + pub fn bezier_elements(&self) -> SplineResult> { + let g = &self.mesh.grid; + let marks = knot_lines(g, &self.blends)?; + let [cs, ct] = frame(self.degree); + let (ns, nt) = (g.ns(), g.nt()); + let (ci_end, cj_end) = (ns - 1 - cs, nt - 1 - ct); + let mut taken = vec![false; ns * nt]; + let is_taken = |taken: &[bool], i: usize, j: usize| taken.get(j * ns + i).copied().unwrap_or(true); + let mut rects = Vec::new(); + for cj in ct..cj_end { + for ci in cs..ci_end { + if is_taken(&taken, ci, cj) { + continue; + } + let mut i1 = ci + 1; + while i1 < ci_end && !is_taken(&taken, i1, cj) && !marks.v(i1, cj) { + i1 += 1; + } + let mut j1 = cj + 1; + while j1 < cj_end && (ci..i1).all(|x| !is_taken(&taken, x, j1) && !marks.h(x, j1)) && (ci + 1..i1).all(|x| !marks.v(x, j1)) { + j1 += 1; + } + for y in cj..j1 { + for x in ci..i1 { + if let Some(c) = taken.get_mut(y * ns + x) { + *c = true; + } + } + } + rects.push((ci, i1, cj, j1)); + } + } + let val = |k: &[f64], i: usize| k.get(i).copied().unwrap_or(0.0); + let [ps, pt] = self.degree; + let mut work = 0usize; + let mut out = Vec::new(); + for (i0, i1, j0, j1) in rects { + let (s0, s1, t0, t1) = (val(&g.s, i0), val(&g.s, i1), val(&g.t, j0), val(&g.t, j1)); + if !(s0 < s1 && t0 < t1) { + continue; + } + let mut ctrl = vec![HPoint::ZERO; (ps + 1) * (pt + 1)]; + for (b, c) in self.blends.iter().zip(&self.ctrl) { + let [[bs0, bs1], [bt0, bt1]] = b.support(); + if !(s1.min(bs1) > s0.max(bs0) && t1.min(bt1) > t0.max(bt0)) { + continue; + } + work += ctrl.len(); + if work > MAX_EXTRACT_WORK { + return Err(SplineError::TooLarge("Bézier extraction needs too much work".into())); + } + let cs = basis_bezier(&b.sv, s0, s1)?; + let ct = basis_bezier(&b.tv, t0, t1)?; + for (j, bj) in ct.iter().enumerate() { + for (i, bi) in cs.iter().enumerate() { + if let Some(slot) = ctrl.get_mut(j * (ps + 1) + i) { + *slot += *c * (bi * bj); + } + } + } + } + if ctrl.iter().any(|h| !h.is_finite() || !(h.w > 0.0)) { + return invalid("a Bézier patch has a non-positive weight"); + } + out.push(BezierPatch { s: [s0, s1], t: [t0, t1], degree: self.degree, ctrl }); + } + Ok(out) + } +} + +/// The unit edges of the index grid that carry a knot line of some blending function: each +/// function's s-knots across its t-support and its t-knots across its s-support. +fn knot_lines(g: &Grid, blends: &[Blend]) -> SplineResult { + let mut m = Grid::empty(g.s.clone(), g.t.clone())?; + let mut work = 0usize; + for b in blends { + let (Some(&s0), Some(&s1), Some(&t0), Some(&t1)) = (b.s.first(), b.s.last(), b.t.first(), b.t.last()) else { continue }; + work += (b.s.len() * (t1 - t0)) + (b.t.len() * (s1 - s0)); + if work > MAX_EXTRACT_WORK { + return Err(SplineError::TooLarge("Bézier extraction needs too much work".into())); + } + for &x in &b.s { + (t0..t1).for_each(|y| m.set_v(x, y)); + } + for &y in &b.t { + (s0..s1).for_each(|x| m.set_h(x, y)); + } + } + Ok(m) +} diff --git a/crates/geom/src/tspline/mesh.rs b/crates/geom/src/tspline/mesh.rs new file mode 100644 index 00000000..85ff816d --- /dev/null +++ b/crates/geom/src/tspline/mesh.rs @@ -0,0 +1,442 @@ +//! The T-mesh in index space (Sederberg, Zheng, Bakenov, Nasri 2003, §3; Scott, Li, +//! Sederberg, Hughes 2012, §2): a rectangle of index lines — columns `0..ns` carrying the +//! non-decreasing knot values `s[i]`, rows `0..nt` carrying `t[j]` — and the unit edges of +//! that grid that are present. Equal knot values on neighbouring lines give zero-length knot +//! intervals (repeated knots), which is how clamped boundaries and multiple knots are written. +//! +//! A valid T-mesh has its outer rectangle, rectangular faces (every grid point is empty, an +//! edge passing straight through, a T-junction, a crossing, or one of the four outer +//! corners), and at least one edge on every index line. A *vertex* is a point on both a column +//! edge and a row edge; every vertex carries a control point. + +use crate::bspline::{Dir, SplineError, SplineResult, check_knots, invalid}; + +/// Most index lines in either direction. +pub const MAX_INDEX_LINES: usize = 1 << 12; +/// Most grid points (`ns · nt`). +pub const MAX_GRID_POINTS: usize = 1 << 20; + +/// A run of edges on one index line, given to [`TMesh::from_segments`]. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Segment { + /// Edges on column `i` (constant s) from row `j0` to row `j1`. + Column { i: usize, j0: usize, j1: usize }, + /// Edges on row `j` (constant t) from column `i0` to column `i1`. + Row { j: usize, i0: usize, i1: usize }, +} + +/// A rectangle in index space, `[i0, i1] × [j0, j1]`. +#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub struct IndexRect { + pub i0: usize, + pub i1: usize, + pub j0: usize, + pub j1: usize, +} + +/// The index lines and the unit edges present on them; no validity implied. +#[derive(Clone, Debug, PartialEq)] +pub(crate) struct Grid { + pub(crate) s: Vec, + pub(crate) t: Vec, + /// Row edge `(i, j)–(i+1, j)` at `j·(ns−1) + i`. + h: Vec, + /// Column edge `(i, j)–(i, j+1)` at `i·(nt−1) + j`. + v: Vec, +} + +impl Grid { + /// An empty grid (no edges) after checking sizes and knots. + pub(crate) fn empty(s: Vec, t: Vec) -> SplineResult { + for k in [&s, &t] { + if k.len() < 2 { + return invalid("a T-mesh needs at least two index lines each way"); + } + if k.len() > MAX_INDEX_LINES { + return Err(SplineError::TooLarge(format!("{} index lines (at most {MAX_INDEX_LINES})", k.len()))); + } + check_knots(k)?; + } + if s.len() * t.len() > MAX_GRID_POINTS { + return Err(SplineError::TooLarge(format!("{} grid points (at most {MAX_GRID_POINTS})", s.len() * t.len()))); + } + let (ns, nt) = (s.len(), t.len()); + Ok(Grid { h: vec![false; (ns - 1) * nt], v: vec![false; ns * (nt - 1)], s, t }) + } + + pub(crate) fn ns(&self) -> usize { + self.s.len() + } + pub(crate) fn nt(&self) -> usize { + self.t.len() + } + + /// Row edge from `(i, j)` to `(i+1, j)`. + pub(crate) fn h(&self, i: usize, j: usize) -> bool { + i + 1 < self.ns() && self.h.get(j * (self.ns() - 1) + i).copied().unwrap_or(false) + } + /// Column edge from `(i, j)` to `(i, j+1)`. + pub(crate) fn v(&self, i: usize, j: usize) -> bool { + j + 1 < self.nt() && self.v.get(i * (self.nt() - 1) + j).copied().unwrap_or(false) + } + pub(crate) fn set_h(&mut self, i: usize, j: usize) { + let ns = self.ns(); + if i + 1 < ns + && let Some(e) = self.h.get_mut(j * (ns - 1) + i) + { + *e = true; + } + } + pub(crate) fn set_v(&mut self, i: usize, j: usize) { + let nt = self.nt(); + if j + 1 < nt + && let Some(e) = self.v.get_mut(i * (nt - 1) + j) + { + *e = true; + } + } + + /// The point `(i, j)` lies on a column edge (passing through or ending there). + pub(crate) fn on_column(&self, i: usize, j: usize) -> bool { + (j > 0 && self.v(i, j - 1)) || self.v(i, j) + } + /// The point `(i, j)` lies on a row edge. + pub(crate) fn on_row(&self, i: usize, j: usize) -> bool { + (i > 0 && self.h(i - 1, j)) || self.h(i, j) + } + + /// The first `count` index lines a ray from `(i, j)` crosses (2003 §3: the knots of a + /// blending function are read off rays cast from its anchor). Along s it crosses columns + /// carrying an edge at row `j`; along t, rows carrying an edge at column `i`. Nearest + /// first; a ray that leaves the mesh repeats the boundary line (clamped ends). + pub(crate) fn ray(&self, i: usize, j: usize, dir: Dir, forward: bool, count: usize) -> Vec { + let (start, last) = dir.pick((i, self.ns() - 1), (j, self.nt() - 1)); + let hit = |x: usize| match dir { + Dir::S => self.on_column(x, j), + Dir::T => self.on_row(i, x), + }; + let mut out = Vec::with_capacity(count); + let mut x = start; + while out.len() < count { + let next = if forward { (x < last).then(|| x + 1) } else { x.checked_sub(1) }; + match next { + Some(n) => { + x = n; + if hit(x) { + out.push(x); + } + } + None => out.push(x), + } + } + out + } + + /// The local knot vector (as line indices) of a degree-`p` blending function anchored at + /// `(i, j)` in direction `dir`: `(p+1)/2` crossings each side of the anchor's own line + /// (odd degree; 2003 §3). + pub(crate) fn local_knots(&self, i: usize, j: usize, dir: Dir, p: usize) -> Vec { + let h = p.div_ceil(2); + let own = dir.pick(i, j); + let mut k: Vec = self.ray(i, j, dir, false, h); + k.reverse(); + k.push(own); + k.extend(self.ray(i, j, dir, true, h)); + k + } + + /// The face containing the unit cell whose lower-left corner is `(ci, cj)`. + pub(crate) fn face_of_cell(&self, ci: usize, cj: usize) -> IndexRect { + let i0 = (0..=ci).rev().find(|&x| self.v(x, cj)).unwrap_or(0); + let i1 = (ci + 1..self.ns()).find(|&x| self.v(x, cj)).unwrap_or(self.ns() - 1); + let j0 = (0..=cj).rev().find(|&y| self.h(ci, y)).unwrap_or(0); + let j1 = (cj + 1..self.nt()).find(|&y| self.h(ci, y)).unwrap_or(self.nt() - 1); + IndexRect { i0, i1, j0, j1 } + } + + /// Split a face with a full-height column edge at `x` so that the point `(x, y)` lies on + /// a column edge (the face above the point, or below it on the top row). + pub(crate) fn split_column_at(&mut self, x: usize, y: usize) { + let cj = if y + 1 < self.nt() { y } else { y.saturating_sub(1) }; + let f = self.face_of_cell(x.min(self.ns() - 2), cj); + for j in f.j0..f.j1 { + self.set_v(x, j); + } + } + /// Split a face with a full-width row edge at `y` so that `(x, y)` lies on a row edge (the + /// face to the right of the point, or to its left on the last column). + pub(crate) fn split_row_at(&mut self, x: usize, y: usize) { + let ci = if x + 1 < self.ns() { x } else { x.saturating_sub(1) }; + let f = self.face_of_cell(ci, y.min(self.nt() - 2)); + for i in f.i0..f.i1 { + self.set_h(i, y); + } + } + + /// Insert a new column at index `at` (`0 < at < ns`) with knot value `value`; row edges + /// crossing it are split in two, the new column has no edges. + pub(crate) fn insert_column(&mut self, at: usize, value: f64) -> SplineResult<()> { + let (ns, nt) = (self.ns(), self.nt()); + if at == 0 || at >= ns || ns + 1 > MAX_INDEX_LINES || (ns + 1) * nt > MAX_GRID_POINTS { + return Err(SplineError::TooLarge("refinement needs too many index lines".into())); + } + let mut h = Vec::with_capacity(ns * nt); + for j in 0..nt { + for i in 0..ns { + // New row edge i runs (i, i+1) in the new numbering. + let old = if i + 1 < at { + self.h(i, j) + } else if i + 1 == at || i == at { + self.h(at - 1, j) + } else { + self.h(i - 1, j) + }; + h.push(old); + } + } + let mut v = Vec::with_capacity((ns + 1) * (nt - 1)); + for i in 0..=ns { + for j in 0..nt - 1 { + v.push(match i.cmp(&at) { + std::cmp::Ordering::Less => self.v(i, j), + std::cmp::Ordering::Equal => false, + std::cmp::Ordering::Greater => self.v(i - 1, j), + }); + } + } + self.s.insert(at, value); + self.h = h; + self.v = v; + Ok(()) + } + + /// The grid transposed (s ↔ t), so row operations can reuse column code. + pub(crate) fn transposed(&self) -> Grid { + let (ns, nt) = (self.ns(), self.nt()); + let mut g = Grid { s: self.t.clone(), t: self.s.clone(), h: vec![false; (nt - 1) * ns], v: vec![false; nt * (ns - 1)] }; + for i in 0..ns { + for j in 0..nt { + if self.h(i, j) { + g.set_v(j, i); + } + if self.v(i, j) { + g.set_h(j, i); + } + } + } + g + } +} + +/// A validated T-mesh. +#[derive(Clone, Debug, PartialEq)] +pub struct TMesh { + pub(crate) grid: Grid, + vertices: Vec<(usize, usize)>, + /// Vertex number of grid point `(i, j)` at `j·ns + i`, `u32::MAX` where there is none. + vid: Vec, +} + +impl TMesh { + /// The full tensor-product mesh on the given knot lines. + pub fn tensor(s: Vec, t: Vec) -> SplineResult { + let (ns, nt) = (s.len(), t.len()); + let segs: Vec = (0..ns) + .map(|i| Segment::Column { i, j0: 0, j1: nt.saturating_sub(1) }) + .chain((0..nt).map(|j| Segment::Row { j, i0: 0, i1: ns.saturating_sub(1) })) + .collect(); + TMesh::from_segments(s, t, &segs) + } + + /// A T-mesh from knot lines and runs of edges (the outer rectangle must be among them). + pub fn from_segments(s: Vec, t: Vec, segments: &[Segment]) -> SplineResult { + let mut g = Grid::empty(s, t)?; + if segments.len() > MAX_GRID_POINTS { + return Err(SplineError::TooLarge("too many segments".into())); + } + let (ns, nt) = (g.ns(), g.nt()); + for seg in segments { + match *seg { + Segment::Column { i, j0, j1 } => { + if i >= ns || j0 >= j1 || j1 >= nt { + return invalid(format!("column segment {i}: rows {j0}..{j1} are not inside the mesh")); + } + (j0..j1).for_each(|j| g.set_v(i, j)); + } + Segment::Row { j, i0, i1 } => { + if j >= nt || i0 >= i1 || i1 >= ns { + return invalid(format!("row segment {j}: columns {i0}..{i1} are not inside the mesh")); + } + (i0..i1).for_each(|i| g.set_h(i, j)); + } + } + } + TMesh::from_grid(g) + } + + /// Validate a grid as a T-mesh. + pub(crate) fn from_grid(g: Grid) -> SplineResult { + let (ns, nt) = (g.ns(), g.nt()); + if (0..ns - 1).any(|i| !g.h(i, 0) || !g.h(i, nt - 1)) || (0..nt - 1).any(|j| !g.v(0, j) || !g.v(ns - 1, j)) { + return invalid("the T-mesh's outer rectangle is not closed"); + } + if let Some(i) = (0..ns).find(|&i| !(0..nt - 1).any(|j| g.v(i, j))) { + return invalid(format!("column {i} carries no edge")); + } + if let Some(j) = (0..nt).find(|&j| !(0..ns - 1).any(|i| g.h(i, j))) { + return invalid(format!("row {j} carries no edge")); + } + let mut vertices = Vec::new(); + let mut vid = vec![u32::MAX; ns * nt]; + for j in 0..nt { + for i in 0..ns { + let l = i > 0 && g.h(i - 1, j); + let r = g.h(i, j); + let d = j > 0 && g.v(i, j - 1); + let u = g.v(i, j); + let n = [l, r, d, u].iter().filter(|&&e| e).count(); + let corner = (i == 0 || i == ns - 1) && (j == 0 || j == nt - 1); + match n { + 1 => return invalid(format!("an edge ends without meeting another at ({i}, {j})")), + 2 if (l != r) && !corner => return invalid(format!("a face is not rectangular at ({i}, {j})")), + _ => {} + } + if (l || r) && (d || u) { + if let Some(slot) = vid.get_mut(j * ns + i) { + *slot = vertices.len() as u32; + } + vertices.push((i, j)); + } + } + } + Ok(TMesh { grid: g, vertices, vid }) + } + + /// Knot values of the columns. + pub fn s_knots(&self) -> &[f64] { + &self.grid.s + } + /// Knot values of the rows. + pub fn t_knots(&self) -> &[f64] { + &self.grid.t + } + /// The vertices `(column, row)`, row by row (column fastest). Control points follow this + /// order. + pub fn vertices(&self) -> &[(usize, usize)] { + &self.vertices + } + /// Vertex number of grid point `(i, j)`. + pub fn vertex_at(&self, i: usize, j: usize) -> Option { + if i >= self.grid.ns() { + return None; + } + self.vid.get(j * self.grid.ns() + i).filter(|&&k| k != u32::MAX).map(|&k| k as usize) + } + /// The parameter point `(s, t)` of vertex `k`. + pub fn vertex_param(&self, k: usize) -> Option<(f64, f64)> { + let &(i, j) = self.vertices.get(k)?; + Some((*self.grid.s.get(i)?, *self.grid.t.get(j)?)) + } + /// Is the edge from `(i, j)` to `(i+1, j)` present? + pub fn row_edge(&self, i: usize, j: usize) -> bool { + self.grid.h(i, j) + } + /// Is the edge from `(i, j)` to `(i, j+1)` present? + pub fn column_edge(&self, i: usize, j: usize) -> bool { + self.grid.v(i, j) + } + /// Every unit edge present: a tensor-product mesh without T-junctions. + pub fn is_regular(&self) -> bool { + (0..self.grid.nt()).all(|j| (0..self.grid.ns() - 1).all(|i| self.grid.h(i, j))) + && (0..self.grid.ns()).all(|i| (0..self.grid.nt() - 1).all(|j| self.grid.v(i, j))) + } + /// The faces, by their lower-left corner. + pub fn faces(&self) -> Vec { + let g = &self.grid; + let mut out = Vec::new(); + for cj in 0..g.nt() - 1 { + for ci in 0..g.ns() - 1 { + if g.v(ci, cj) && g.h(ci, cj) { + out.push(g.face_of_cell(ci, cj)); + } + } + } + out + } + + /// The local knot vector (line indices) of a degree-`p` blending function anchored at + /// vertex `(i, j)` in direction `dir` (2003 §3). + pub fn local_knots(&self, i: usize, j: usize, dir: Dir, p: usize) -> Vec { + self.grid.local_knots(i, j, dir, p) + } + + /// The T-junction extensions for bi-degree `degree` (Scott, Li, Sederberg, Hughes 2012, + /// §2, after Li, Zheng, Sederberg, Hughes, Scott 2012): from each interior T-junction a face + /// extension runs in the direction of the missing edge across `⌈p/2⌉` perpendicular lines + /// carrying an edge there, and an edge extension the other way across `⌊p/2⌋` of them. + /// `p` is the degree along the extension. + pub fn extensions(&self, degree: [usize; 2]) -> Vec { + let g = &self.grid; + let (ns, nt) = (g.ns(), g.nt()); + let mut out = Vec::new(); + for &(i, j) in &self.vertices { + if i == 0 || j == 0 || i == ns - 1 || j == nt - 1 { + continue; + } + let (l, r, d, u) = (g.h(i - 1, j), g.h(i, j), g.v(i, j - 1), g.v(i, j)); + if [l, r, d, u].iter().filter(|&&e| e).count() != 3 { + continue; + } + let (dir, missing_forward) = if !l || !r { (Dir::S, !r) } else { (Dir::T, !u) }; + let p = dir.pick(degree[0], degree[1]); + let reach = |forward: bool, count: usize| -> usize { + if count == 0 { + return dir.pick(i, j); + } + g.ray(i, j, dir, forward, count).last().copied().unwrap_or(0) + }; + let face = reach(missing_forward, p.div_ceil(2)); + let edge = reach(!missing_forward, p / 2); + let range = if face < edge { [face, edge] } else { [edge, face] }; + let at = dir.pick(j, i); + out.push(Extension { dir, at, range, from: (i, j) }); + } + out + } + + /// Analysis-suitable for bi-degree `degree`: no extension along s meets an extension along + /// t (Scott et al. 2012, §2). Analysis-suitable T-splines have linearly independent + /// blending functions that form a partition of unity. + pub fn is_analysis_suitable(&self, degree: [usize; 2]) -> bool { + self.conflict(degree).is_none() + } + + /// The first extension along s that meets an extension along t. + pub(crate) fn conflict(&self, degree: [usize; 2]) -> Option { + let ns = self.grid.ns(); + let ext = self.extensions(degree); + let mut covered: Vec> = vec![None; ns * self.grid.nt()]; + for (k, e) in ext.iter().enumerate().filter(|(_, e)| e.dir == Dir::S) { + for x in e.range[0]..=e.range[1] { + if let Some(c) = covered.get_mut(e.at * ns + x) { + *c = Some(k); + } + } + } + ext.iter() + .filter(|e| e.dir == Dir::T) + .find_map(|e| (e.range[0]..=e.range[1]).find_map(|y| covered.get(y * ns + e.at).copied().flatten())) + .and_then(|k| ext.get(k).copied()) + } +} + +/// A T-junction extension: a closed segment along `dir` on index line `at` (the row for +/// [`Dir::S`], the column for [`Dir::T`]) spanning the lines `range`, from the T-junction +/// `from`. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub struct Extension { + pub dir: Dir, + pub at: usize, + pub range: [usize; 2], + pub from: (usize, usize), +} diff --git a/crates/geom/src/tspline/mod.rs b/crates/geom/src/tspline/mod.rs new file mode 100644 index 00000000..531f13b5 --- /dev/null +++ b/crates/geom/src/tspline/mod.rs @@ -0,0 +1,219 @@ +//! T-spline surfaces: control points on a T-mesh, which allows T-junctions, so a surface can +//! be refined locally and patches with different knot vectors join without gaps. +//! +//! References: Sederberg, Zheng, Bakenov, Nasri 2003, "T-splines and T-NURCCs" (ACM ToG +//! 22(3)) — §3 the T-spline equation and the local knot vectors read off the T-mesh by rays; +//! Sederberg, Cardon, Finnigan, North, Zheng, Lyche 2004, "T-spline simplification and local +//! refinement" (ACM ToG 23(3)) — §3 local refinement by knot insertion into blending +//! functions; Scott, Li, Sederberg, Hughes 2012, "Local refinement of analysis-suitable +//! T-splines" (CMAME 213–216) — §2 index-space T-meshes and T-junction extensions; Piegl & +//! Tiller, *The NURBS Book* — the univariate basis. +//! +//! Scope: odd degree in each direction (control points at vertices, as in the 2003 paper; +//! even degrees put anchors on faces and are refused), rectangular index domains. +//! +//! The surface is `S(s, t) = Σ wₖ Pₖ Bₖ(s, t) / Σ wₖ Bₖ(s, t)` with `Bₖ(s, t) = N[sₖ](s) · +//! N[tₖ](t)`, the univariate B-splines of the vertex's local knot vectors (2003, eq. 2–3). +//! Near the boundary the outermost `(p−1)/2` index lines each way form a frame that must be +//! complete lines; the domain is the rectangle inside it. + +mod extract; +mod mesh; +mod refine; + +pub use crate::bspline::Dir; +pub use mesh::{Extension, IndexRect, MAX_GRID_POINTS, MAX_INDEX_LINES, Segment, TMesh}; +pub use refine::Cut; + +use crate::Vec3; +use crate::bspline::{HPoint, MAX_SPLINE_DEGREE, Side, SplineError, SplineResult, invalid, local_basis, weighted_points}; + +/// A blending function: its anchor's local knot vectors as line indices and as values, and +/// its support. +#[derive(Clone, Debug, PartialEq)] +struct Blend { + s: Vec, + t: Vec, + sv: Vec, + tv: Vec, +} + +impl Blend { + fn support(&self) -> [[f64; 2]; 2] { + let ends = |v: &[f64]| [v.first().copied().unwrap_or(0.0), v.last().copied().unwrap_or(0.0)]; + [ends(&self.sv), ends(&self.tv)] + } +} + +/// The number of frame lines outside the domain on each side, `(p−1)/2` per direction. +fn frame(degree: [usize; 2]) -> [usize; 2] { + degree.map(|p| p.saturating_sub(1) / 2) +} + +/// A T-spline surface of odd bi-degree. +#[derive(Clone, Debug, PartialEq)] +pub struct TSpline { + mesh: TMesh, + degree: [usize; 2], + ctrl: Vec, + blends: Vec, +} + +impl TSpline { + /// A T-spline from a mesh, a bi-degree and one point and positive weight per vertex (in + /// [`TMesh::vertices`] order). + pub fn new(mesh: TMesh, degree: [usize; 2], points: &[Vec3], weights: &[f64]) -> SplineResult { + TSpline::from_homogeneous(mesh, degree, weighted_points(points, weights)?) + } + + /// A T-spline from homogeneous control points (`(x·w, y·w, z·w, w)`, `w > 0`). + pub fn from_homogeneous(mesh: TMesh, degree: [usize; 2], ctrl: Vec) -> SplineResult { + for p in degree { + if p == 0 || p > MAX_SPLINE_DEGREE { + return Err(SplineError::Unsupported(format!("degree {p} (1 to {MAX_SPLINE_DEGREE})"))); + } + if p % 2 == 0 { + return Err(SplineError::Unsupported(format!("even degree {p} (T-splines here have odd degrees)"))); + } + } + let g = &mesh.grid; + let (ns, nt) = (g.ns(), g.nt()); + let [cs, ct] = frame(degree); + if ns < 2 * cs + 2 || nt < 2 * ct + 2 { + return invalid("too few index lines for the degree"); + } + let full_col = |i: usize| (0..nt - 1).all(|j| g.v(i, j)); + let full_row = |j: usize| (0..ns - 1).all(|i| g.h(i, j)); + if (0..=cs).any(|k| !full_col(k) || !full_col(ns - 1 - k)) || (0..=ct).any(|k| !full_row(k) || !full_row(nt - 1 - k)) { + return invalid("the outermost (p−1)/2 index lines must be complete"); + } + if !(g.s.get(cs) < g.s.get(ns - 1 - cs)) || !(g.t.get(ct) < g.t.get(nt - 1 - ct)) { + return invalid("the T-spline's domain is empty"); + } + if ctrl.len() != mesh.vertices().len() { + return invalid(format!("{} control points for {} vertices", ctrl.len(), mesh.vertices().len())); + } + if ctrl.iter().any(|h| !h.is_finite() || !(h.w > 0.0) || h.project().is_none()) { + return invalid("control points must be finite with positive weights"); + } + let mut blends = Vec::with_capacity(ctrl.len()); + for &(i, j) in mesh.vertices() { + let s = g.local_knots(i, j, Dir::S, degree[0]); + let t = g.local_knots(i, j, Dir::T, degree[1]); + let sv: Vec = s.iter().map(|&x| g.s.get(x).copied().unwrap_or(0.0)).collect(); + let tv: Vec = t.iter().map(|&y| g.t.get(y).copied().unwrap_or(0.0)).collect(); + let b = Blend { s, t, sv, tv }; + let [[s0, s1], [t0, t1]] = b.support(); + if !(s0 < s1 && t0 < t1) { + return invalid(format!("repeated knots leave the blending function at ({i}, {j}) empty")); + } + blends.push(b); + } + Ok(TSpline { mesh, degree, ctrl, blends }) + } + + pub fn mesh(&self) -> &TMesh { + &self.mesh + } + pub fn degree(&self) -> [usize; 2] { + self.degree + } + /// Homogeneous control points, one per vertex. + pub fn homogeneous(&self) -> &[HPoint] { + &self.ctrl + } + /// Cartesian control points, one per vertex. + pub fn points(&self) -> Vec { + self.ctrl.iter().map(|h| h.project().unwrap_or(Vec3::ZERO)).collect() + } + /// Weights, one per vertex. + pub fn weights(&self) -> Vec { + self.ctrl.iter().map(|h| h.w).collect() + } + /// The local knot vectors (values) of vertex `k`'s blending function. + pub fn local_knots(&self, k: usize) -> Option<(&[f64], &[f64])> { + self.blends.get(k).map(|b| (b.sv.as_slice(), b.tv.as_slice())) + } + + /// The parameter domain `[[s₀, s₁], [t₀, t₁]]`. + pub fn domain(&self) -> [[f64; 2]; 2] { + let g = &self.mesh.grid; + let [cs, ct] = frame(self.degree); + let at = |k: &[f64], i: usize| k.get(i).copied().unwrap_or(0.0); + [[at(&g.s, cs), at(&g.s, g.ns() - 1 - cs)], [at(&g.t, ct), at(&g.t, g.nt() - 1 - ct)]] + } + + /// Analysis-suitable (no s-extension meets a t-extension; Scott et al. 2012, §2). + pub fn is_analysis_suitable(&self) -> bool { + self.mesh.is_analysis_suitable(self.degree) + } + + /// Which side of the knots to evaluate on, after checking `(s, t)` lies in the domain. + fn sides(&self, s: f64, t: f64) -> SplineResult<[Side; 2]> { + let [ds, dt] = self.domain(); + if !(s >= ds[0] && s <= ds[1] && t >= dt[0] && t <= dt[1]) { + return invalid(format!("({s}, {t}) is outside the T-spline's domain")); + } + let side = |x: f64, hi: f64| if x >= hi { Side::Left } else { Side::Right }; + Ok([side(s, ds[1]), side(t, dt[1])]) + } + + /// Blending function `k` at `(s, t)`, with its first derivatives `[B, ∂B/∂s, ∂B/∂t]`. + fn blend_at(&self, b: &Blend, s: f64, t: f64, sides: [Side; 2], n: usize) -> [f64; 3] { + let [[s0, s1], [t0, t1]] = b.support(); + if s < s0 || s > s1 || t < t0 || t > t1 { + return [0.0; 3]; + } + let ns = local_basis(&b.sv, s, sides[0], n); + let nt = local_basis(&b.tv, t, sides[1], n); + let (n0, n1) = (ns.first().copied().unwrap_or(0.0), ns.get(1).copied().unwrap_or(0.0)); + let (m0, m1) = (nt.first().copied().unwrap_or(0.0), nt.get(1).copied().unwrap_or(0.0)); + [n0 * m0, n1 * m0, n0 * m1] + } + + /// The value of vertex `k`'s blending function `Bₖ(s, t)` (no weight). + pub fn blending(&self, k: usize, s: f64, t: f64) -> SplineResult { + let sides = self.sides(s, t)?; + let b = self.blends.get(k).ok_or_else(|| SplineError::Invalid(format!("no vertex {k}")))?; + Ok(self.blend_at(b, s, t, sides, 0)[0]) + } + + /// Homogeneous sums `[A, ∂A/∂s, ∂A/∂t]` at `(s, t)`. + fn sums(&self, s: f64, t: f64, n: usize) -> SplineResult<[HPoint; 3]> { + let sides = self.sides(s, t)?; + let mut acc = [HPoint::ZERO; 3]; + for (b, c) in self.blends.iter().zip(&self.ctrl) { + let v = self.blend_at(b, s, t, sides, n); + for (a, x) in acc.iter_mut().zip(v) { + if x != 0.0 { + *a += *c * x; + } + } + } + Ok(acc) + } + + /// The point at `(s, t)`. + pub fn eval(&self, s: f64, t: f64) -> SplineResult { + let [a, _, _] = self.sums(s, t, 0)?; + a.project().ok_or_else(|| SplineError::Invalid(format!("the blending functions vanish at ({s}, {t})"))) + } + + /// The point and the first partial derivatives `[S, ∂S/∂s, ∂S/∂t]` at `(s, t)` (quotient + /// rule on the homogeneous sums). + pub fn derivs(&self, s: f64, t: f64) -> SplineResult<[Vec3; 3]> { + let [a, a_s, a_t] = self.sums(s, t, 1)?; + let p = a.project().ok_or_else(|| SplineError::Invalid(format!("the blending functions vanish at ({s}, {t})")))?; + let d = |h: HPoint| (h.xyz() - p * h.w) / a.w; + Ok([p, d(a_s), d(a_t)]) + } + + /// The unit normal `∂S/∂s × ∂S/∂t` at `(s, t)`, if the surface is regular there. + pub fn normal(&self, s: f64, t: f64) -> SplineResult> { + let [_, ds, dt] = self.derivs(s, t)?; + Ok(ds.cross(dt).normalized()) + } +} + +#[cfg(test)] +mod tests; diff --git a/crates/geom/src/tspline/refine.rs b/crates/geom/src/tspline/refine.rs new file mode 100644 index 00000000..75f3af85 --- /dev/null +++ b/crates/geom/src/tspline/refine.rs @@ -0,0 +1,317 @@ +//! Local refinement that keeps the surface exactly (Sederberg, Cardon, Finnigan, North, +//! Zheng, Lyche 2004, §3). +//! +//! The requested knot lines are cut into the T-mesh first. Then every blending function of +//! the original surface is carried as a record — its local knot vectors (as index lines) and +//! its homogeneous coefficient — and the records are made to agree with the refined mesh: +//! +//! * **Violation 1** — the mesh's ray from a record's anchor crosses a line inside the +//! record's support that the record lacks: insert that knot into the record, which splits +//! it into two records (knot insertion into one blending function, 2004 §3). +//! * **Violation 2** — the record has a knot where the ray crosses nothing: split the face +//! there so the mesh dictates that knot. +//! * **Violation 3** — a record's anchor is not a vertex: split a face so it becomes one. +//! +//! Records with the same knot vectors merge. When nothing is violated, each record is exactly +//! the blending function of the vertex at its anchor, so summing the coefficients per vertex +//! gives control points of the refined T-spline, and the surface is unchanged. All knots are +//! handled as line indices, so every comparison is exact; no new index line appears after +//! the cuts, so the work is bounded by the grid. + +use std::collections::BTreeMap; + +use super::TSpline; +use super::mesh::{Grid, TMesh}; +use crate::bspline::{Dir, HPoint, SplineError, SplineResult, invalid, split_basis, split_knots}; + +/// Most cuts in one refinement. +const MAX_CUTS: usize = 1 << 12; +/// Most record checks and mesh changes in one refinement. +const MAX_REFINE_STEPS: usize = 1 << 22; + +/// A knot line to cut into a T-spline. +#[derive(Clone, Copy, Debug, PartialEq)] +pub enum Cut { + /// The line `s = s` from `t = t0` to `t = t1`: every face it crosses there is split. + S { s: f64, t0: f64, t1: f64 }, + /// The line `t = t` from `s = s0` to `s = s1`. + T { t: f64, s0: f64, s1: f64 }, +} + +impl Cut { + /// The direction the cut line is constant in: [`Dir::S`] for a line `s = …`. + fn dir(&self) -> Dir { + match self { + Cut::S { .. } => Dir::S, + Cut::T { .. } => Dir::T, + } + } +} + +/// A record's knot vectors (s, t) as line indices. +type Key = (Vec, Vec); + +enum Fix { + /// Violation 1: insert line `k` into the record along the direction. + Split(Dir, usize), + /// Violations 2 and 3: the point (column, row) must lie on a column edge. + Column(usize, usize), + /// … on a row edge. + Row(usize, usize), +} + +impl TSpline { + /// The same surface on a refined T-mesh with the given knot lines cut in (2004 §3; see the + /// module notes). The cut lines must lie strictly inside the domain; a cut at the value + /// of an existing knot line extends that line. Further faces may be split to keep the + /// surface exact and, when this T-spline is analysis-suitable, to keep the result + /// analysis-suitable (T-junctions are carried across faces until no s-extension meets a + /// t-extension; Scott et al. 2012). + pub fn refine(&self, cuts: &[Cut]) -> SplineResult { + if cuts.len() > MAX_CUTS { + return Err(SplineError::TooLarge(format!("{} cuts (at most {MAX_CUTS})", cuts.len()))); + } + let mut g = self.mesh.grid.clone(); + let mut smap: Vec = (0..g.ns()).collect(); + let mut tmap: Vec = (0..g.nt()).collect(); + let [ds, dt] = self.domain(); + // Row cuts run as column cuts on the transposed grid; consecutive cuts along one + // direction share a transposition. + for run in cuts.chunk_by(|a, b| a.dir() == b.dir()) { + let along_t = run.first().is_some_and(|c| c.dir() == Dir::T); + if along_t { + g = g.transposed(); + } + for cut in run { + match *cut { + Cut::S { s, t0, t1 } => cut_column(&mut g, &mut smap, s, [t0, t1], ds)?, + Cut::T { t, s0, s1 } => cut_column(&mut g, &mut tmap, t, [s0, s1], dt)?, + } + } + if along_t { + g = g.transposed(); + } + } + let remap = |idx: &[usize], map: &[usize]| -> Vec { idx.iter().map(|&x| map.get(x).copied().unwrap_or(x)).collect() }; + let mut rec: BTreeMap = BTreeMap::new(); + for (b, c) in self.blends.iter().zip(&self.ctrl) { + *rec.entry((remap(&b.s, &smap), remap(&b.t, &tmap))).or_insert(HPoint::ZERO) += *c; + } + // Keep analysis-suitability (Scott et al. 2012): while an s-extension meets a + // t-extension, carry the s-extension's T-junction across its face, then make the + // records agree with the mesh again. Every round adds edges, so this ends. + let keep_suitable = self.is_analysis_suitable(); + let mesh = loop { + close(&mut g, &mut rec, self.degree)?; + let mesh = TMesh::from_grid(g.clone())?; + match mesh.conflict(self.degree).filter(|_| keep_suitable) { + Some(e) => extend_t_junction(&mut g, e.from), + None => break mesh, + } + }; + let mut ctrl = vec![HPoint::ZERO; mesh.vertices().len()]; + let mut seen = vec![false; ctrl.len()]; + for ((s, t), c) in &rec { + let (ci, cj) = (centre(s), centre(t)); + let k = mesh.vertex_at(ci, cj).ok_or_else(|| SplineError::Invalid("refinement left a blending function without a vertex".into()))?; + if let (Some(slot), Some(f)) = (ctrl.get_mut(k), seen.get_mut(k)) { + *slot += *c; + *f = true; + } + } + if let Some(k) = seen.iter().position(|&f| !f) { + let (i, j) = mesh.vertices().get(k).copied().unwrap_or_default(); + return Err(SplineError::Unsupported(format!( + "this refinement leaves the control point at ({i}, {j}) without a blending function (the refined T-mesh is not analysis-suitable)" + ))); + } + let out = TSpline::from_homogeneous(mesh, self.degree, ctrl)?; + // The records are the refined blending functions, one per vertex. + for (s, t) in rec.keys() { + let k = out.mesh.vertex_at(centre(s), centre(t)); + if k.and_then(|k| out.blends.get(k)).is_none_or(|b| &b.s != s || &b.t != t) { + return Err(SplineError::Invalid("refinement did not converge to the refined mesh's blending functions".into())); + } + } + Ok(out) + } +} + +/// Carry a T-junction's missing edge across the face beyond it. +fn extend_t_junction(g: &mut Grid, (i, j): (usize, usize)) { + let (l, r, d, u) = (i > 0 && g.h(i - 1, j), g.h(i, j), j > 0 && g.v(i, j - 1), g.v(i, j)); + if !r { + g.split_row_at(i, j); + } else if !l { + g.split_row_at(i.saturating_sub(1), j); + } else if !u { + g.split_column_at(i, j); + } else if !d { + g.split_column_at(i, j.saturating_sub(1)); + } +} + +/// The anchor line of a local knot vector. +fn centre(k: &[usize]) -> usize { + k.get(k.len() / 2).copied().unwrap_or(0) +} + +/// Cut the column line `value` into every face it crosses between `range[0]` and `range[1]`, +/// adding the column if needed. +fn cut_column(g: &mut Grid, map: &mut [usize], value: f64, range: [f64; 2], dom: [f64; 2]) -> SplineResult<()> { + let [a, b] = range; + if !(value.is_finite() && a.is_finite() && b.is_finite()) { + return invalid("a cut needs finite numbers"); + } + if !(a < b) { + return invalid("a cut needs a non-empty range"); + } + if !(value > dom[0] && value < dom[1]) { + return invalid(format!("the cut at {value} is not inside the domain")); + } + let same: Vec = (0..g.ns()).filter(|&x| g.s.get(x) == Some(&value)).collect(); + let (x, new) = match same.as_slice() { + [] => { + let at = g.s.iter().position(|&k| k > value).unwrap_or(g.ns() - 1); + g.insert_column(at, value)?; + map.iter_mut().filter(|m| **m >= at).for_each(|m| *m += 1); + (at, true) + } + [x] => (*x, false), + _ => return Err(SplineError::Unsupported(format!("cutting along the repeated knot {value}"))), + }; + let mut added = false; + for y in 0..g.nt() - 1 { + if g.v(x, y) { + continue; + } + let f = g.face_of_cell(x, y); + let (fa, fb) = (g.t.get(f.j0).copied().unwrap_or(0.0), g.t.get(f.j1).copied().unwrap_or(0.0)); + // Faces overlapping the range, and faces of zero height inside it (so a cut from a + // clamped domain end runs on through the zero-height frame to the boundary). + if fb.min(b) - fa.max(a) > 0.0 || (fa == fb && fa >= a && fb <= b) { + (f.j0..f.j1).for_each(|j| g.set_v(x, j)); + added = true; + } + } + if new && !added { + return invalid(format!("the cut at {value} crosses no face")); + } + Ok(()) +} + +/// Resolve the violations until every record is a blending function of the mesh. +fn close(g: &mut Grid, rec: &mut BTreeMap, degree: [usize; 2]) -> SplineResult<()> { + let mut steps = 0usize; + let mut tick = |n: usize| -> SplineResult<()> { + steps += n; + if steps > MAX_REFINE_STEPS { Err(SplineError::TooLarge("the refinement needs too much work".into())) } else { Ok(()) } + }; + loop { + // Violation 1 first (knot insertion into records), until none is left. + loop { + let keys: Vec = rec.keys().cloned().collect(); + tick(keys.len())?; + let mut changed = false; + for key in keys { + if let Some(Fix::Split(dir, k)) = check(g, &key, degree) + && let Some(c) = rec.remove(&key) + { + for (piece, f) in split(g, &key, dir, k)? { + if f != 0.0 { + *rec.entry(piece).or_insert(HPoint::ZERO) += c * f; + } + } + changed = true; + } + } + if !changed { + break; + } + } + // Then one mesh change for the first record that needs one. + tick(rec.len())?; + let fix = rec.keys().find_map(|k| match check(g, k, degree) { + Some(Fix::Split(..)) | None => None, + other => other, + }); + match fix { + Some(Fix::Column(x, y)) => g.split_column_at(x, y), + Some(Fix::Row(x, y)) => g.split_row_at(x, y), + _ => return Ok(()), + } + } +} + +/// The first violation of a record against the mesh, if any. +fn check(g: &Grid, key: &Key, degree: [usize; 2]) -> Option { + let (s, t) = key; + let (ci, cj) = (centre(s), centre(t)); + let ms = g.local_knots(ci, cj, Dir::S, degree[0]); + match compare(s, &ms) { + Some(Ok(k)) => return Some(Fix::Split(Dir::S, k)), + Some(Err(x)) => return Some(Fix::Column(x, cj)), + None => {} + } + let mt = g.local_knots(ci, cj, Dir::T, degree[1]); + match compare(t, &mt) { + Some(Ok(k)) => return Some(Fix::Split(Dir::T, k)), + Some(Err(y)) => return Some(Fix::Row(ci, y)), + None => {} + } + if !g.on_column(ci, cj) { + return Some(Fix::Column(ci, cj)); + } + if !g.on_row(ci, cj) { + return Some(Fix::Row(ci, cj)); + } + None +} + +/// Compare a record's knots with the mesh's, outward from the anchor: `Ok(k)` the mesh +/// crosses line `k` first (insert it into the record), `Err(x)` the record has line `x` where +/// the mesh crosses nothing (the mesh must gain it). +fn compare(rec: &[usize], mesh: &[usize]) -> Option> { + let mid = rec.len() / 2; + for k in mid + 1..rec.len() { + let (r, m) = (rec.get(k)?, mesh.get(k)?); + if m < r { + return Some(Ok(*m)); + } + if m > r { + return Some(Err(*r)); + } + } + for k in (0..mid).rev() { + let (r, m) = (rec.get(k)?, mesh.get(k)?); + if m > r { + return Some(Ok(*m)); + } + if m < r { + return Some(Err(*r)); + } + } + None +} + +/// Insert line `k` into a record along `dir`: the two pieces and their factors (2004 §3). +/// Pieces whose knots all have one value are the zero function and are dropped. +fn split(g: &Grid, key: &Key, dir: Dir, k: usize) -> SplineResult> { + let (idx, vals) = dir.pick((&key.0, &g.s), (&key.1, &g.t)); + let value = |x: usize| vals.get(x).copied().unwrap_or(0.0); + let values: Vec = idx.iter().map(|&x| value(x)).collect(); + let (c1, c2) = split_basis(&values, value(k)).ok_or_else(|| SplineError::Invalid("refinement inserted a knot outside a support".into()))?; + let mut out = Vec::with_capacity(2); + for (piece, c) in split_knots(idx, k).into_iter().zip([c1, c2]) { + let (Some(&lo), Some(&hi)) = (piece.first(), piece.last()) else { continue }; + if value(lo) == value(hi) { + continue; + } + let key = match dir { + Dir::S => (piece, key.1.clone()), + Dir::T => (key.0.clone(), piece), + }; + out.push((key, c)); + } + Ok(out) +} diff --git a/crates/geom/src/tspline/tests.rs b/crates/geom/src/tspline/tests.rs new file mode 100644 index 00000000..9791d78a --- /dev/null +++ b/crates/geom/src/tspline/tests.rs @@ -0,0 +1,350 @@ +use super::*; +use crate::bspline::tests::{Lcg, random_surface}; +use crate::bspline::{BezierPatch, NurbsSurface, SplineError}; + +fn uniform(n: usize) -> Vec { + (0..n).map(|k| k as f64).collect() +} + +fn col(i: usize, j0: usize, j1: usize) -> Segment { + Segment::Column { i, j0, j1 } +} +fn row(j: usize, i0: usize, i1: usize) -> Segment { + Segment::Row { j, i0, i1 } +} + +/// 7 × 7 lines, row 3 only from column 0 to column 3: one T-junction at (3, 3) whose +/// extension runs along s. +fn t_mesh() -> TMesh { + let mut segs: Vec = (0..7).map(|i| col(i, 0, 6)).collect(); + segs.extend([0, 1, 2, 4, 5, 6].map(|j| row(j, 0, 6))); + segs.push(row(3, 0, 3)); + TMesh::from_segments(uniform(7), uniform(7), &segs).unwrap() +} + +fn random_points(seed: u64, n: usize) -> (Vec, Vec) { + let mut rng = Lcg(seed); + let pts = (0..n).map(|_| Vec3::new(rng.unit(), rng.unit(), rng.unit()) * 10.0).collect(); + let w = (0..n).map(|_| 0.5 + rng.unit()).collect(); + (pts, w) +} + +fn t_spline(seed: u64) -> TSpline { + let mesh = t_mesh(); + let (pts, w) = random_points(seed, mesh.vertices().len()); + TSpline::new(mesh, [3, 3], &pts, &w).unwrap() +} + +/// A clamped bicubic NURBS with uneven interior knots, as a T-spline. +fn tensor_spline(seed: u64) -> (NurbsSurface, TSpline) { + let ku = vec![0.0, 0.0, 0.0, 0.0, 0.2, 0.45, 0.5, 0.8, 1.0, 1.0, 1.0, 1.0]; + let kv = vec![0.0, 0.0, 0.0, 0.0, 0.3, 0.6, 1.0, 1.0, 1.0, 1.0]; + let n = random_surface(seed, ku, kv, [3, 3]); + let ts = TSpline::from_nurbs(&n).unwrap(); + (n, ts) +} + +fn grid_params(dom: [[f64; 2]; 2], k: usize) -> Vec<(f64, f64)> { + let mut out = Vec::new(); + for i in 0..=k { + for j in 0..=k { + let s = dom[0][0] + (dom[0][1] - dom[0][0]) * i as f64 / k as f64; + let t = dom[1][0] + (dom[1][1] - dom[1][0]) * j as f64 / k as f64; + out.push((s, t)); + } + } + out +} + +fn max_deviation(a: &TSpline, b: &TSpline) -> f64 { + assert_eq!(a.domain(), b.domain()); + grid_params(a.domain(), 60).into_iter().map(|(s, t)| a.eval(s, t).unwrap().dist(b.eval(s, t).unwrap())).fold(0.0, f64::max) +} + +#[test] +fn tmesh_validation_refuses_malformed_meshes() { + let frame = |extra: &[Segment]| -> Vec { + let mut s = vec![col(0, 0, 4), col(4, 0, 4), row(0, 0, 4), row(4, 0, 4)]; + s.extend_from_slice(extra); + s + }; + let full: Vec = (1..4).flat_map(|k| [col(k, 0, 4), row(k, 0, 4)]).collect(); + assert!(TMesh::from_segments(uniform(5), uniform(5), &frame(&full)).is_ok()); + // A dangling edge, an L-shaped corner, an open boundary, an empty line. + let mut dangling = full.clone(); + dangling.push(col(2, 0, 4)); + dangling.retain(|s| *s != row(2, 0, 4)); + dangling.push(row(2, 0, 1)); + dangling.push(row(2, 3, 4)); + assert!(TMesh::from_segments(uniform(5), uniform(5), &frame(&dangling)).is_ok(), "T-junctions on both sides are fine"); + let dangling_edge = frame(&[col(1, 0, 4), col(3, 0, 4), row(1, 0, 4), row(2, 0, 4), row(3, 0, 1), col(2, 0, 3)]); + assert!(TMesh::from_segments(uniform(5), uniform(5), &dangling_edge).is_err(), "edge ending inside a face"); + let mut l_shape = frame(&[col(1, 0, 4), col(2, 0, 4), col(3, 0, 4), row(1, 0, 4), row(3, 0, 4)]); + l_shape.push(row(2, 0, 2)); + l_shape.retain(|s| *s != col(2, 0, 4)); + l_shape.extend([col(2, 0, 2), col(2, 3, 4)]); + assert!(TMesh::from_segments(uniform(5), uniform(5), &l_shape).is_err(), "L-shaped face corner"); + assert!(TMesh::from_segments(uniform(5), uniform(5), &full).is_err(), "no outer rectangle"); + let mut empty = frame(&full); + empty.retain(|s| *s != col(2, 0, 4)); + assert!(TMesh::from_segments(uniform(5), uniform(5), &empty).is_err(), "a column without edges"); + assert!(TMesh::from_segments(vec![0.0, 1.0, f64::NAN, 3.0, 4.0], uniform(5), &frame(&full)).is_err()); + assert!(TMesh::from_segments(vec![0.0, 2.0, 1.0, 3.0, 4.0], uniform(5), &frame(&full)).is_err()); + assert!(TMesh::from_segments(uniform(5), uniform(5), &[col(9, 0, 4)]).is_err()); + assert!(TMesh::from_segments(uniform(5), uniform(5), &[col(1, 3, 3)]).is_err()); + assert!(matches!(TMesh::tensor(uniform(MAX_INDEX_LINES + 1), uniform(2)), Err(SplineError::TooLarge(_)))); + assert!(matches!(TMesh::tensor(uniform(2048), uniform(2048)), Err(SplineError::TooLarge(_)))); + assert!(TMesh::tensor(uniform(1), uniform(3)).is_err()); +} + +#[test] +fn local_knot_vectors_are_read_off_rays() { + let ts = t_spline(1); + let m = ts.mesh(); + assert!(m.vertex_at(4, 3).is_none(), "(4, 3) lies on a column edge only"); + assert!(m.vertex_at(3, 3).is_some(), "the T-junction is a vertex"); + let knots = |i, j| { + let k = m.vertex_at(i, j).unwrap(); + let (s, t) = ts.local_knots(k).unwrap(); + (s.to_vec(), t.to_vec()) + }; + assert_eq!(knots(2, 3), (vec![0.0, 1.0, 2.0, 3.0, 4.0], vec![1.0, 2.0, 3.0, 4.0, 5.0])); + assert_eq!(knots(3, 4).1, vec![2.0, 3.0, 4.0, 5.0, 6.0]); + // The ray down column 4 passes (4, 3), which carries no row edge. + assert_eq!(knots(4, 4).1, vec![1.0, 2.0, 4.0, 5.0, 6.0]); + assert_eq!(knots(4, 2).1, vec![0.0, 1.0, 2.0, 4.0, 5.0]); + // Rays leaving the mesh repeat the boundary line. + assert_eq!(knots(0, 0), (vec![0.0, 0.0, 0.0, 1.0, 2.0], vec![0.0, 0.0, 0.0, 1.0, 2.0])); + assert_eq!(knots(6, 5).1, vec![2.0, 4.0, 5.0, 6.0, 6.0]); +} + +#[test] +fn regular_tspline_is_the_tensor_nurbs_surface() { + for (seed, degree, ku, kv) in [ + (2, [3, 3], vec![0.0, 0.0, 0.0, 0.0, 0.2, 0.45, 0.5, 0.8, 1.0, 1.0, 1.0, 1.0], vec![0.0, 0.0, 0.0, 0.0, 0.3, 0.6, 1.0, 1.0, 1.0, 1.0]), + (3, [3, 1], vec![0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 1.0, 1.0, 1.0, 1.0], vec![0.0, 0.0, 0.4, 1.0, 1.0]), + (4, [5, 3], vec![0.0; 6].into_iter().chain([0.3, 0.7]).chain([1.0; 6]).collect(), vec![0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0, 1.0]), + ] { + let n = random_surface(seed, ku, kv, degree); + let ts = TSpline::from_nurbs(&n).unwrap(); + assert!(ts.mesh().is_regular() && ts.is_analysis_suitable()); + assert_eq!(ts.domain(), n.domain()); + let worst = grid_params(n.domain(), 40).into_iter().map(|(s, t)| ts.eval(s, t).unwrap().dist(n.eval(s, t).unwrap())).fold(0.0, f64::max); + assert!(worst < 1e-12, "degree {degree:?}: T-spline differs from NURBS by {worst}"); + assert_eq!(ts.to_nurbs().unwrap(), n, "NURBS → T-spline → NURBS is exact"); + } + // An unclamped NURBS is clamped on the way in (same surface). + let n = random_surface(5, uniform(10), uniform(9), [3, 3]); + let ts = TSpline::from_nurbs(&n).unwrap(); + let worst = grid_params(n.domain(), 30).into_iter().map(|(s, t)| ts.eval(s, t).unwrap().dist(n.eval(s, t).unwrap())).fold(0.0, f64::max); + assert!(worst < 1e-12, "unclamped: {worst}"); + assert!(TSpline::from_nurbs(&random_surface(6, uniform(9), uniform(9), [2, 3])).is_err(), "even degree"); +} + +#[test] +fn blending_functions_sum_to_one_with_t_junctions() { + let ts = t_spline(1); + assert!(ts.is_analysis_suitable()); + let n = ts.mesh().vertices().len(); + let mut worst = 0.0f64; + for (s, t) in grid_params(ts.domain(), 64) { + let sum: f64 = (0..n).map(|k| ts.blending(k, s, t).unwrap()).sum(); + worst = worst.max((sum - 1.0).abs()); + } + assert!(worst < 1e-14, "partition of unity off by {worst}"); +} + +#[test] +fn derivatives_match_finite_differences() { + let ts = t_spline(9); + let h = 1e-6; + let mut worst = 0.0f64; + for (s, t) in grid_params([[1.05, 4.95], [1.05, 4.95]], 17) { + let [_, ds, dt] = ts.derivs(s, t).unwrap(); + let fs = (ts.eval(s + h, t).unwrap() - ts.eval(s - h, t).unwrap()) / (2.0 * h); + let ft = (ts.eval(s, t + h).unwrap() - ts.eval(s, t - h).unwrap()) / (2.0 * h); + worst = worst.max((ds - fs).len()).max((dt - ft).len()); + } + assert!(worst < 1e-6, "derivative error {worst}"); + assert!(ts.normal(3.0, 3.0).unwrap().is_some()); +} + +#[test] +fn analysis_suitability_follows_t_junction_extensions() { + assert!(tensor_spline(1).1.is_analysis_suitable()); + let ext = t_mesh().extensions([3, 3]); + assert_eq!(ext, vec![Extension { dir: Dir::S, at: 3, range: [2, 5], from: (3, 3) }]); + assert!(t_mesh().is_analysis_suitable([3, 3])); + // Row 4 stops at column 4 (extension along s over columns 3..7 on row 4: column 5 has no + // edge there) and column 5 stops at row 3 (extension along t over rows 2..6 on column 5): + // they cross at (5, 4). + let mut segs: Vec = [0, 1, 2, 3, 4, 6, 7, 8].iter().map(|&i| col(i, 0, 8)).collect(); + segs.push(col(5, 0, 3)); + segs.extend([0, 1, 2, 3, 5, 6, 7, 8].map(|j| row(j, 0, 8))); + segs.push(row(4, 0, 4)); + let m = TMesh::from_segments(uniform(9), uniform(9), &segs).unwrap(); + let ext = m.extensions([3, 3]); + assert!(ext.contains(&Extension { dir: Dir::S, at: 4, range: [3, 7], from: (4, 4) })); + assert!(ext.contains(&Extension { dir: Dir::T, at: 5, range: [2, 6], from: (5, 3) })); + assert!(!m.is_analysis_suitable([3, 3])); + // Refining it does not repair it (only suitable T-splines are kept suitable), but stays exact. + let (pts, w) = random_points(4, m.vertices().len()); + let ts = TSpline::new(m, [3, 3], &pts, &w).unwrap(); + let r = ts.refine(&[Cut::S { s: 2.5, t0: 2.0, t1: 3.0 }]).unwrap(); + assert!(max_deviation(&ts, &r) < 1e-12); +} + +#[test] +fn a_cut_across_the_domain_is_boehm_knot_insertion() { + let (n, ts) = tensor_spline(4); + let r = ts.refine(&[Cut::S { s: 0.33, t0: 0.0, t1: 1.0 }, Cut::T { t: 0.75, s0: 0.0, s1: 1.0 }]).unwrap(); + assert!(r.mesh().is_regular()); + let want = n.insert_knot(Dir::S, 0.33).unwrap().insert_knot(Dir::T, 0.75).unwrap(); + let got = r.to_nurbs().unwrap(); + assert_eq!(got.knots(), want.knots()); + let worst = got.homogeneous().iter().zip(want.homogeneous()).map(|(a, b)| (*a - *b).xyz().len() + (a.w - b.w).abs()).fold(0.0, f64::max); + assert!(worst < 1e-12, "control points differ by {worst}"); +} + +#[test] +fn local_refinement_keeps_the_surface() { + let (_, ts) = tensor_spline(7); + let before = ts.mesh().vertices().len(); + let r = ts.refine(&[Cut::S { s: 0.33, t0: 0.3, t1: 0.6 }]).unwrap(); + let dev = max_deviation(&ts, &r); + eprintln!("local cut: {} → {} vertices, AS {}, max deviation {dev:e}", before, r.mesh().vertices().len(), r.is_analysis_suitable()); + assert!(dev < 1e-12, "deviation {dev}"); + assert!(r.mesh().vertices().len() <= before + 4, "a one-face cut stays local"); + assert!(r.is_analysis_suitable()); + // Refining again near the new T-junctions, and across them. + let r2 = r.refine(&[Cut::T { t: 0.45, s0: 0.2, s1: 0.45 }, Cut::S { s: 0.39, t0: 0.3, t1: 0.45 }]).unwrap(); + let dev2 = max_deviation(&ts, &r2); + eprintln!("second refinement: {} vertices, max deviation {dev2:e}", r2.mesh().vertices().len()); + assert!(dev2 < 1e-12, "deviation {dev2}"); + assert!(r2.is_analysis_suitable(), "refinement keeps analysis-suitability"); + let n2 = r2.mesh().vertices().len(); + let pu = grid_params(r2.domain(), 50) + .into_iter() + .map(|(s, t)| ((0..n2).map(|k| r2.blending(k, s, t).unwrap()).sum::() - 1.0).abs()) + .fold(0.0, f64::max); + assert!(pu < 1e-14, "partition of unity after refinement off by {pu}"); + // A rational T-spline with T-junctions. + let t = t_spline(3); + let r3 = t.refine(&[Cut::S { s: 2.5, t0: 2.0, t1: 4.0 }, Cut::T { t: 3.5, s0: 3.0, s1: 5.0 }]).unwrap(); + let dev3 = max_deviation(&t, &r3); + eprintln!("rational T-spline: {} → {} vertices, max deviation {dev3:e}", t.mesh().vertices().len(), r3.mesh().vertices().len()); + assert!(dev3 < 1e-12, "deviation {dev3}"); +} + +#[test] +fn random_local_refinements_keep_the_surface() { + let mut rng = Lcg(42); + let mut worst = 0.0f64; + let mut vertices = 0; + for round in 0..40 { + let (_, ts) = tensor_spline(round); + let mut cur = ts.clone(); + for _ in 0..1 + (rng.unit() * 3.0) as usize { + let v = 0.05 + 0.9 * (rng.unit() * 40.0).round() / 40.0; + let a = (rng.unit() * 8.0).round() / 10.0; + let b = (a + 0.1 + (rng.unit() * 5.0).round() / 10.0).min(1.0); + let cut = if rng.unit() < 0.5 { Cut::S { s: v, t0: a, t1: b } } else { Cut::T { t: v, s0: a, s1: b } }; + cur = cur.refine(&[cut]).unwrap_or_else(|e| panic!("round {round}: {cut:?}: {e}")); + } + assert!(cur.is_analysis_suitable(), "round {round}: refinement keeps analysis-suitability"); + vertices += cur.mesh().vertices().len() - ts.mesh().vertices().len(); + worst = worst.max(max_deviation(&ts, &cur)); + } + eprintln!("random refinements: max deviation {worst:e}, {vertices} vertices added over 40 rounds"); + assert!(worst < 1e-12, "deviation {worst}"); +} + +fn check_patches(ts: &TSpline, patches: &[BezierPatch]) -> f64 { + let [ds, dt] = ts.domain(); + let area: f64 = patches.iter().map(|p| (p.s[1] - p.s[0]) * (p.t[1] - p.t[0])).sum(); + assert!((area - (ds[1] - ds[0]) * (dt[1] - dt[0])).abs() < 1e-12, "the patches tile the domain"); + let mut worst = 0.0f64; + for p in patches { + for (s, t) in grid_params([p.s, p.t], 6) { + worst = worst.max(p.eval(s, t).unwrap().dist(ts.eval(s, t).unwrap())); + } + } + worst +} + +#[test] +fn bezier_extraction_reproduces_the_surface() { + let ts = t_spline(5); + let patches = ts.bezier_elements().unwrap(); + let dev = check_patches(&ts, &patches); + eprintln!("Bézier extraction, T-junction mesh: {} patches, max deviation {dev:e}", patches.len()); + assert!(dev < 1e-12); + let (n, tensor) = tensor_spline(8); + let patches = tensor.bezier_elements().unwrap(); + assert_eq!(patches.len(), 5 * 3, "one patch per pair of knot spans"); + assert!(check_patches(&tensor, &patches) < 1e-12); + let _ = n; + let refined = tensor.refine(&[Cut::S { s: 0.33, t0: 0.3, t1: 0.6 }, Cut::T { t: 0.45, s0: 0.2, s1: 0.45 }]).unwrap(); + let patches = refined.bezier_elements().unwrap(); + let dev = check_patches(&refined, &patches); + eprintln!("Bézier extraction, refined: {} patches, max deviation {dev:e}", patches.len()); + assert!(dev < 1e-12); +} + +#[test] +fn hostile_tspline_input_is_refused() { + let mesh = t_mesh(); + let n = mesh.vertices().len(); + let (pts, w) = random_points(1, n); + assert!(TSpline::new(mesh.clone(), [3, 3], &pts, &w).is_ok()); + let mut bad = pts.clone(); + bad[3] = Vec3::new(f64::NAN, 0.0, 0.0); + assert!(TSpline::new(mesh.clone(), [3, 3], &bad, &w).is_err()); + let mut bw = w.clone(); + bw[0] = 0.0; + assert!(TSpline::new(mesh.clone(), [3, 3], &pts, &bw).is_err()); + bw[0] = f64::INFINITY; + assert!(TSpline::new(mesh.clone(), [3, 3], &pts, &bw).is_err()); + assert!(TSpline::new(mesh.clone(), [3, 3], &pts[1..], &w[1..]).is_err()); + assert!(matches!(TSpline::new(mesh.clone(), [2, 3], &pts, &w), Err(SplineError::Unsupported(_)))); + assert!(TSpline::new(mesh.clone(), [0, 3], &pts, &w).is_err()); + assert!(TSpline::new(mesh.clone(), [9, 3], &pts, &w).is_err()); + // Degree 5 needs the three outermost lines each way complete; row 3 is not, so a mesh + // whose row 2 stops short is refused. + let mut segs: Vec = (0..7).map(|i| col(i, 0, 6)).collect(); + segs.extend([0, 1, 3, 4, 5, 6].map(|j| row(j, 0, 6))); + segs.push(row(2, 0, 3)); + let short_frame = TMesh::from_segments(uniform(7), uniform(7), &segs).unwrap(); + let (p6, w6) = random_points(3, short_frame.vertices().len()); + assert!(TSpline::new(short_frame.clone(), [3, 3], &p6, &w6).is_ok()); + assert!(TSpline::new(short_frame, [5, 5], &p6, &w6).is_err()); + // A row of knots repeated so often that a blending function is empty. + let m = TMesh::tensor(vec![0.0, 0.0, 0.0, 1.0, 2.0], uniform(5)).unwrap(); + let (p, w5) = random_points(2, 25); + assert!(TSpline::new(m, [3, 3], &p, &w5).is_err()); + // An empty domain. + let m = TMesh::tensor(vec![0.0, 1.0, 1.0, 2.0], uniform(4)).unwrap(); + assert!(TSpline::new(m, [3, 3], &p[..16], &w5[..16]).is_err()); + + let ts = TSpline::new(mesh, [3, 3], &pts, &w).unwrap(); + for (s, t) in [(f64::NAN, 2.0), (2.0, f64::INFINITY), (0.5, 2.0), (2.0, 5.5)] { + assert!(ts.eval(s, t).is_err() && ts.derivs(s, t).is_err() && ts.blending(0, s, t).is_err()); + } + assert!(ts.blending(n + 5, 2.0, 2.0).is_err()); + for cut in [ + Cut::S { s: f64::NAN, t0: 1.0, t1: 2.0 }, + Cut::S { s: 2.5, t0: 2.0, t1: 1.0 }, + Cut::S { s: 2.5, t0: 2.0, t1: f64::INFINITY }, + Cut::S { s: 0.5, t0: 1.0, t1: 2.0 }, + Cut::T { t: 5.0, s0: 1.0, s1: 2.0 }, + Cut::T { t: 9.0, s0: 1.0, s1: 2.0 }, + Cut::S { s: 2.5, t0: 8.0, t1: 9.0 }, + ] { + assert!(ts.refine(&[cut]).is_err(), "{cut:?}"); + } + assert!(matches!(ts.refine(&vec![Cut::S { s: 2.5, t0: 1.0, t1: 2.0 }; 5000]), Err(SplineError::TooLarge(_)))); + // Refining the same line twice is a no-op the second time. + let a = ts.refine(&[Cut::S { s: 2.5, t0: 1.0, t1: 5.0 }]).unwrap(); + let b = a.refine(&[Cut::S { s: 2.5, t0: 1.0, t1: 5.0 }]).unwrap(); + assert_eq!(a, b); +} diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index 2d914411..a0defffa 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -64,6 +64,7 @@ mod sweep_path; #[cfg(test)] mod thread_tests; mod topo; +mod tspline; pub use blend::{ChamferSide, chamfer, chamfer_sides, fillet}; pub use blend_corner::fillet_radii; @@ -96,6 +97,7 @@ pub use surfaces::{copy_faces, extend, patch_edges, patch_region, stitch, thicke pub use sweep_path::sweep_path; pub use topo::{CylinderFace, cylinder_face_at}; pub use topo::{TopoCounts, merged_topology, seam_flags}; +pub use tspline::tspline_sheet; #[derive(Debug, Clone, PartialEq, thiserror::Error)] pub enum KernelError { diff --git a/crates/kernel/src/surfaces.rs b/crates/kernel/src/surfaces.rs index 3101dccc..360f23ac 100644 --- a/crates/kernel/src/surfaces.rs +++ b/crates/kernel/src/surfaces.rs @@ -22,7 +22,7 @@ impl Body { } /// A body of faces as they are (no orientation fix: a surface has no inside). -fn surface_body(faces: Vec) -> Body { +pub(crate) fn surface_body(faces: Vec) -> Body { surface_of(mt::Shell::from(faces)) } diff --git a/crates/kernel/src/tspline.rs b/crates/kernel/src/tspline.rs new file mode 100644 index 00000000..e804765c --- /dev/null +++ b/crates/kernel/src/tspline.rs @@ -0,0 +1,267 @@ +//! T-spline sheets: a T-spline surface (`solvecraft_geom::tspline`) as a surface body of NURBS +//! faces. A T-mesh without T-junctions becomes one clamped tensor face; otherwise each Bézier +//! element is a face, and the faces share their edges and vertices — a side that meets two +//! neighbours is split at the T-point between them — so the sheet is one sewn open shell. + +use std::collections::hash_map::Entry; +use std::collections::{BTreeMap, BTreeSet, HashMap}; + +use solvecraft_geom::bspline::{BezierPatch, HPoint, NurbsSurface, SplineError, bezier_segment}; +use solvecraft_geom::tspline::TSpline; +use truck_modeling as mt; + +use crate::body::Body; +use crate::{KernelError, Result, guard}; + +/// A T-spline as a surface body (an open shell of NURBS faces whose normals follow +/// `∂S/∂s × ∂S/∂t`). +pub fn tspline_sheet(ts: &TSpline) -> Result { + let faces = if ts.mesh().is_regular() { + let n = ts.to_nurbs().and_then(|n| n.clamped()).map_err(spline_err)?; + guard("t-spline sheet", || tensor_face(&n).map(|f| vec![f]))? + } else { + let patches = ts.bezier_elements().map_err(spline_err)?; + guard("t-spline sheet", || patch_faces(&patches))? + }; + Ok(crate::surfaces::surface_body(faces)) +} + +fn spline_err(e: SplineError) -> KernelError { + KernelError::Invalid(format!("t-spline: {e}")) +} + +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(); + let [ku, kv] = n.knots(); + let (ku, kv) = (mt::KnotVec::from(ku.clone()), mt::KnotVec::from(kv.clone())); + let corners = [n.point(0, 0), n.point(nu - 1, 0), n.point(nu - 1, nv - 1), n.point(0, nv - 1)]; + let v: Vec = corners.into_iter().map(vertex).collect::>()?; + let row = |j: usize| -> Vec { (0..nu).map(|i| n.point(i, j)).collect() }; + let column = |i: usize| -> Vec { (0..nv).map(|j| n.point(i, j)).collect() }; + let (Some(v0), Some(v1), Some(v2), Some(v3)) = (v.first(), v.get(1), v.get(2), v.get(3)) else { return Err(fail("corners")) }; + let edge = |a: &mt::Vertex, b: &mt::Vertex, c: mt::Curve| mt::Edge::try_new(a, b, c).map_err(fail); + let bottom = edge(v0, v1, nurbs_curve(ku.clone(), &row(0))?)?; + let right = edge(v1, v2, nurbs_curve(kv.clone(), &column(nu - 1))?)?; + 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) +} + +/// A parameter point, exact (knot values are shared, so equal points have equal bits). +type Key = (u64, u64); + +fn key(s: f64, t: f64) -> Key { + ((s + 0.0).to_bits(), (t + 0.0).to_bits()) +} + +/// One face per Bézier patch, sharing vertices and edges; sides are split where a +/// neighbour's corner lands on them. +fn patch_faces(patches: &[BezierPatch]) -> Result> { + // Corner parameters on each line of constant t (their s values) and of constant s. + let mut on_t: BTreeMap> = BTreeMap::new(); + let mut on_s: BTreeMap> = BTreeMap::new(); + for p in patches { + for s in p.s { + for t in p.t { + let (ks, kt) = key(s, t); + on_t.entry(kt).or_default().insert(ks); + on_s.entry(ks).or_default().insert(kt); + } + } + } + let mut verts: HashMap = HashMap::new(); + let mut edges: HashMap<(Key, Key), mt::Edge> = HashMap::new(); + let mut faces = Vec::with_capacity(patches.len()); + for p in patches { + let mut wire: Vec = Vec::new(); + for side in 0..4 { + // Side 0: t = t₀, 1: s = s₁, 2: t = t₁, 3: s = s₀; each runs with its parameter + // rising, sides 2 and 3 are walked backwards. + let (along_s, fixed, [lo, hi]) = match side { + 0 => (true, p.t[0], p.s), + 1 => (false, p.s[1], p.t), + 2 => (true, p.t[1], p.s), + _ => (false, p.s[0], p.t), + }; + let line = if along_s { on_t.get(&key(0.0, fixed).1) } else { on_s.get(&key(fixed, 0.0).0) }; + let mut breaks: Vec = line.into_iter().flatten().map(|&b| f64::from_bits(b)).filter(|&x| x > lo && x < hi).collect(); + breaks.insert(0, lo); + breaks.push(hi); + let ctrl = p.side(side); + let at = |x: f64| if along_s { key(x, fixed) } else { key(fixed, x) }; + let mut pieces = Vec::with_capacity(breaks.len()); + for w in breaks.windows(2) { + let &[x0, x1] = w else { continue }; + let (ka, kb) = (at(x0), at(x1)); + let e = if let Some(e) = edges.get(&(ka, kb)) { + e.clone() + } else if let Some(e) = edges.get(&(kb, ka)) { + e.inverse() + } else { + let seg = bezier_segment(&ctrl, (x0 - lo) / (hi - lo), (x1 - lo) / (hi - lo)).map_err(spline_err)?; + let (Some(&first), Some(&last)) = (seg.first(), seg.last()) else { return Err(fail("an empty side")) }; + let mut corner = |k: Key, h: HPoint| -> Result { + Ok(match verts.entry(k) { + Entry::Occupied(e) => e.get().clone(), + Entry::Vacant(e) => e.insert(vertex(h)?).clone(), + }) + }; + let (va, vb) = (corner(ka, first)?, corner(kb, last)?); + let curve = nurbs_curve(mt::KnotVec::bezier_knot(seg.len() - 1), &seg)?; + let e = mt::Edge::try_new(&va, &vb, curve).map_err(fail)?; + edges.insert((ka, kb), e.clone()); + e + }; + pieces.push(e); + } + if side >= 2 { + pieces.reverse(); + pieces.iter_mut().for_each(|e| *e = e.inverse()); + } + wire.extend(pieces); + } + 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)?; + faces.push(mt::Face::try_new(vec![wire.into()], surf).map_err(fail)?); + } + Ok(faces) +} + +#[cfg(test)] +mod tests { + use solvecraft_geom::Vec3; + use solvecraft_geom::bspline::NurbsSurface; + use solvecraft_geom::tspline::{Cut, TSpline}; + use truck_modeling as mt; + + use super::tspline_sheet; + use crate::fuzz::Rng; + + /// A rational bicubic with uneven knots over a 10 × 8 mm area, bulging up to ~3 mm. + fn surface(seed: u64) -> TSpline { + let mut rng = Rng::new(seed); + let ku = vec![0.0, 0.0, 0.0, 0.0, 0.2, 0.45, 0.5, 0.8, 1.0, 1.0, 1.0, 1.0]; + let kv = vec![0.0, 0.0, 0.0, 0.0, 0.3, 0.6, 1.0, 1.0, 1.0, 1.0]; + let (nu, nv) = (ku.len() - 4, kv.len() - 4); + let mut pts = Vec::new(); + for j in 0..nv { + for i in 0..nu { + pts.push(Vec3::new(10.0 * i as f64 / (nu - 1) as f64, 8.0 * j as f64 / (nv - 1) as f64, 3.0 * rng.unit())); + } + } + let w: Vec = (0..pts.len()).map(|_| 0.7 + 0.6 * rng.unit()).collect(); + TSpline::from_nurbs(&NurbsSurface::new([3, 3], [ku, kv], &pts, &w).unwrap()).unwrap() + } + + /// Distance from `p` to the T-spline: nearest of a sample grid, then Gauss–Newton. + fn distance(ts: &TSpline, p: Vec3) -> f64 { + let [ds, dt] = ts.domain(); + let n = 24; + let mut best = (f64::MAX, 0.0, 0.0); + for i in 0..=n { + for j in 0..=n { + let (s, t) = (ds[0] + (ds[1] - ds[0]) * i as f64 / n as f64, dt[0] + (dt[1] - dt[0]) * j as f64 / n as f64); + let d = ts.eval(s, t).unwrap().dist(p); + if d < best.0 { + best = (d, s, t); + } + } + } + let (_, mut s, mut t) = best; + for _ in 0..60 { + let [q, a, b] = ts.derivs(s, t).unwrap(); + 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-30 { + break; + } + let (ga, gb) = (a.dot(r), b.dot(r)); + s = (s - (bb * ga - ab * gb) / det).clamp(ds[0], ds[1]); + t = (t - (aa * gb - ab * ga) / det).clamp(dt[0], dt[1]); + } + ts.eval(s, t).unwrap().dist(p) + } + + fn check_sheet(ts: &TSpline, faces: usize) { + let body = tspline_sheet(ts).unwrap(); + assert!(body.is_surface()); + assert_eq!(body.shell_count(), 1); + assert_eq!(body.face_count(), faces); + let shell = &body.solid.boundaries()[0]; + assert_eq!(shell.shell_condition(), mt::ShellCondition::Oriented, "neighbours share edges, run opposite ways"); + let euler = body.vertex_count() as i64 - body.edge_count() as i64 + body.face_count() as i64; + assert_eq!(euler, 1, "one sewn sheet (a disc)"); + let mesh = body.tessellate(0.01).unwrap(); + assert!(!mesh.triangles.is_empty()); + // Mesh points lie on the surface (not merely within the tessellation tolerance). + let step = (mesh.positions.len() / 300).max(1); + let worst = mesh.positions.iter().step_by(step).map(|p| distance(ts, *p)).fold(0.0, f64::max); + eprintln!("t-spline sheet: {faces} faces, {} mesh points, farthest checked {worst:e} mm from the surface", mesh.positions.len()); + assert!(worst < 1e-6, "a mesh point is {worst} mm off the T-spline"); + } + + #[test] + fn regular_tspline_becomes_one_nurbs_face() { + check_sheet(&surface(1), 1); + } + + #[test] + fn tspline_with_t_junctions_becomes_sewn_bezier_faces() { + let ts = surface(2).refine(&[Cut::S { s: 0.33, t0: 0.3, t1: 0.6 }, Cut::T { t: 0.45, s0: 0.2, s1: 0.45 }]).unwrap(); + assert!(!ts.mesh().is_regular()); + let n = ts.bezier_elements().unwrap().len(); + check_sheet(&ts, n); + } + + #[test] + fn hostile_tsplines_never_panic() { + let base = surface(3); + let n = base.mesh().vertices().len(); + let mesh = base.mesh().clone(); + // Every control point at one place, enormous coordinates, a needle-thin patch. + let cases: Vec> = vec![ + vec![Vec3::new(1.0, 2.0, 3.0); n], + (0..n).map(|k| Vec3::new(1e300, -1e300, k as f64 * 1e299)).collect(), + (0..n).map(|k| Vec3::new(k as f64 * 1e-12, 0.0, 0.0)).collect(), + ]; + for pts in cases { + let Ok(ts) = TSpline::new(mesh.clone(), [3, 3], &pts, &vec![1.0; n]) else { continue }; + if let Ok(body) = tspline_sheet(&ts) { + let _ = body.tessellate(0.01); + } + if let Ok(r) = ts.refine(&[Cut::S { s: 0.33, t0: 0.3, t1: 0.6 }]) + && let Ok(body) = tspline_sheet(&r) + { + let _ = body.tessellate(0.01); + } + } + } +}