* 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 <[email protected]> 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 <[email protected]> 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 <[email protected]> Claude-Session: https://claude.ai/code/session_01J53sR6gEaVt56soFDpTck8 --------- Co-authored-by: Claude Opus 5.5 <[email protected]>
This commit is contained in:
Iraklis KourtisandClaude Opus 5.5 authored and GitHub committed
1 parent
8e87a37511
commit
a04ef4ca59
12 files changed
+2566
-3
No files matched your search
Generated
+1
@@ -3351,6 +3351,7 @@ version = "0.4.0"
|
||||
dependencies = [
|
||||
"serde",
|
||||
"serde_json",
|
||||
"thiserror 2.0.21",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
|
||||
@@ -10,6 +10,7 @@ workspace = true
|
||||
|
||||
[dependencies]
|
||||
serde = { workspace = true }
|
||||
thiserror = { workspace = true }
|
||||
|
||||
[dev-dependencies]
|
||||
serde_json = { workspace = true }
|
||||
@@ -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<T> = Result<T, SplineError>;
|
||||
|
||||
pub(crate) fn invalid<T>(msg: impl Into<String>) -> SplineResult<T> {
|
||||
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<Vec3> {
|
||||
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<f64> 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<X>(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<Vec<HPoint>> {
|
||||
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<f64> {
|
||||
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<T: Copy + PartialOrd>(knots: &[T], k: T) -> [Vec<T>; 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<Vec<f64>> {
|
||||
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<usize> {
|
||||
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<f64>, Vec<HPoint>)> {
|
||||
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<Vec<HPoint>> {
|
||||
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<HPoint>, Vec<HPoint>) {
|
||||
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<HPoint>,
|
||||
}
|
||||
|
||||
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<Vec3> {
|
||||
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<HPoint> {
|
||||
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<f64>; 2],
|
||||
count: [usize; 2],
|
||||
ctrl: Vec<HPoint>,
|
||||
}
|
||||
|
||||
impl NurbsSurface {
|
||||
/// A surface from Cartesian points and positive weights (one per point, u fastest).
|
||||
pub fn new(degree: [usize; 2], knots: [Vec<f64>; 2], points: &[Vec3], weights: &[f64]) -> SplineResult<NurbsSurface> {
|
||||
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<f64>; 2], ctrl: Vec<HPoint>) -> SplineResult<NurbsSurface> {
|
||||
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<f64>; 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<Vec3> {
|
||||
let dom = self.domain();
|
||||
let mut spans = [0; 2];
|
||||
let mut basis: [Vec<f64>; 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<NurbsSurface> {
|
||||
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<HPoint>> = Vec::with_capacity(lines);
|
||||
for line in 0..lines {
|
||||
let row: Vec<HPoint> = (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<NurbsSurface> {
|
||||
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<NurbsSurface> {
|
||||
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<f64> = 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<f64> {
|
||||
(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<f64> = (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<HPoint> = (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<f64>, kv: Vec<f64>, degree: [usize; 2]) -> NurbsSurface {
|
||||
let mut rng = Lcg(seed);
|
||||
let n = (ku.len() - degree[0] - 1) * (kv.len() - degree[1] - 1);
|
||||
let pts: Vec<Vec3> = (0..n).map(|_| Vec3::new(rng.unit(), rng.unit(), rng.unit()) * 10.0).collect();
|
||||
let w: Vec<f64> = (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<f64> = (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<HPoint> = (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]);
|
||||
}
|
||||
}
|
||||
@@ -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};
|
||||
|
||||
@@ -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<TSpline> {
|
||||
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<f64> {
|
||||
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<NurbsSurface> {
|
||||
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<f64> {
|
||||
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<Vec<BezierPatch>> {
|
||||
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<Grid> {
|
||||
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)
|
||||
}
|
||||
@@ -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<f64>,
|
||||
pub(crate) t: Vec<f64>,
|
||||
/// Row edge `(i, j)–(i+1, j)` at `j·(ns−1) + i`.
|
||||
h: Vec<bool>,
|
||||
/// Column edge `(i, j)–(i, j+1)` at `i·(nt−1) + j`.
|
||||
v: Vec<bool>,
|
||||
}
|
||||
|
||||
impl Grid {
|
||||
/// An empty grid (no edges) after checking sizes and knots.
|
||||
pub(crate) fn empty(s: Vec<f64>, t: Vec<f64>) -> SplineResult<Grid> {
|
||||
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<usize> {
|
||||
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<usize> {
|
||||
let h = p.div_ceil(2);
|
||||
let own = dir.pick(i, j);
|
||||
let mut k: Vec<usize> = 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<u32>,
|
||||
}
|
||||
|
||||
impl TMesh {
|
||||
/// The full tensor-product mesh on the given knot lines.
|
||||
pub fn tensor(s: Vec<f64>, t: Vec<f64>) -> SplineResult<TMesh> {
|
||||
let (ns, nt) = (s.len(), t.len());
|
||||
let segs: Vec<Segment> = (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<f64>, t: Vec<f64>, segments: &[Segment]) -> SplineResult<TMesh> {
|
||||
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<TMesh> {
|
||||
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<usize> {
|
||||
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<IndexRect> {
|
||||
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<usize> {
|
||||
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<Extension> {
|
||||
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<Extension> {
|
||||
let ns = self.grid.ns();
|
||||
let ext = self.extensions(degree);
|
||||
let mut covered: Vec<Option<usize>> = 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),
|
||||
}
|
||||
@@ -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<usize>,
|
||||
t: Vec<usize>,
|
||||
sv: Vec<f64>,
|
||||
tv: Vec<f64>,
|
||||
}
|
||||
|
||||
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<HPoint>,
|
||||
blends: Vec<Blend>,
|
||||
}
|
||||
|
||||
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> {
|
||||
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<HPoint>) -> SplineResult<TSpline> {
|
||||
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<f64> = s.iter().map(|&x| g.s.get(x).copied().unwrap_or(0.0)).collect();
|
||||
let tv: Vec<f64> = 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<Vec3> {
|
||||
self.ctrl.iter().map(|h| h.project().unwrap_or(Vec3::ZERO)).collect()
|
||||
}
|
||||
/// Weights, one per vertex.
|
||||
pub fn weights(&self) -> Vec<f64> {
|
||||
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<f64> {
|
||||
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<Vec3> {
|
||||
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<Option<Vec3>> {
|
||||
let [_, ds, dt] = self.derivs(s, t)?;
|
||||
Ok(ds.cross(dt).normalized())
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests;
|
||||
@@ -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<usize>, Vec<usize>);
|
||||
|
||||
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<TSpline> {
|
||||
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<usize> = (0..g.ns()).collect();
|
||||
let mut tmap: Vec<usize> = (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<usize> { idx.iter().map(|&x| map.get(x).copied().unwrap_or(x)).collect() };
|
||||
let mut rec: BTreeMap<Key, HPoint> = 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<usize> = (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<Key, HPoint>, 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<Key> = 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<Fix> {
|
||||
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<Result<usize, usize>> {
|
||||
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<Vec<(Key, f64)>> {
|
||||
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<f64> = 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)
|
||||
}
|
||||
@@ -0,0 +1,350 @@
|
||||
use super::*;
|
||||
use crate::bspline::tests::{Lcg, random_surface};
|
||||
use crate::bspline::{BezierPatch, NurbsSurface, SplineError};
|
||||
|
||||
fn uniform(n: usize) -> Vec<f64> {
|
||||
(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<Segment> = (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<Vec3>, Vec<f64>) {
|
||||
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<Segment> {
|
||||
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<Segment> = (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<Segment> = [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::<f64>() - 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<Segment> = (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);
|
||||
}
|
||||
@@ -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 {
|
||||
|
||||
@@ -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<mt::Face>) -> Body {
|
||||
pub(crate) fn surface_body(faces: Vec<mt::Face>) -> Body {
|
||||
surface_of(mt::Shell::from(faces))
|
||||
}
|
||||
|
||||
|
||||
@@ -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<Body> {
|
||||
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<mt::Vertex> {
|
||||
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<mt::Curve> {
|
||||
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<Vec<mt::Vector4>>) -> Result<mt::Surface> {
|
||||
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<mt::Face> {
|
||||
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<mt::Vertex> = corners.into_iter().map(vertex).collect::<Result<_>>()?;
|
||||
let row = |j: usize| -> Vec<HPoint> { (0..nu).map(|i| n.point(i, j)).collect() };
|
||||
let column = |i: usize| -> Vec<HPoint> { (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<Vec<mt::Vector4>> = (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<Vec<mt::Face>> {
|
||||
// Corner parameters on each line of constant t (their s values) and of constant s.
|
||||
let mut on_t: BTreeMap<u64, BTreeSet<u64>> = BTreeMap::new();
|
||||
let mut on_s: BTreeMap<u64, BTreeSet<u64>> = 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<Key, mt::Vertex> = 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<mt::Edge> = 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<f64> = 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<mt::Vertex> {
|
||||
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<Vec<mt::Vector4>> = (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<f64> = (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<Vec3>> = 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);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
Reference in new issue
Block a user