mirror of
https://github.com/bherbruck/solidcraft.git
synced 2026-11-03 02:25:55 +00:00
geom, kernel: the patch fill meets its neighbours with G0 to G4, on rational edges too
One order-n implementation in the linkage-curve formulation (Pegna & Wolter 1992; Hermann, Lukacs & Wolter 1999): geom::linkage reads the neighbour as a height over its tangent plane to order n (its jet, from the surface's exact partials of any order) and gives each order's condition on the linkage curve, with exact gradients (dual numbers), so fill_iso solves orders 2..n by Newton's method. Each order's row is weighed by h^(k-1) and its jump is held, carried into the fill by Taylor's theorem, within the angle tolerance. The iso fill takes rational edges (a circle's arc on its NURBS form): homogeneous rows, each edge scaled so the corner weights agree, the inside weights the Coons blend and fixed, so the fill stays a linear least-squares problem. Its degree is at least n + 1, refinement adds simple knots, breaks of different sides within rounding are one, and the weak Hermite rows act in the tangent plane at the Coons start's own speed, so they no longer pull against the G1 rows. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LUXeKMoxpM7Vuv8Rb6E3K2
This commit is contained in:
jyviko
@@ -91,6 +91,10 @@ pub struct Across {
|
||||
pub normal: Vec3,
|
||||
/// G2: the neighbour's curvature vector along `leave` (unit speed).
|
||||
pub bend: Option<Vec3>,
|
||||
/// G2 and up: the neighbour as a height over its tangent plane to the order asked
|
||||
/// ([`crate::linkage::Jet`], on the frame `normal × leave`, `leave`), which the conditions of
|
||||
/// [`crate::fill_iso`] read. `None` where only `bend` is known.
|
||||
pub jet: Option<crate::linkage::Jet>,
|
||||
}
|
||||
|
||||
/// One sample of the loop.
|
||||
@@ -226,7 +230,9 @@ pub fn fill_general(
|
||||
fn finite_station(s: &Station) -> bool {
|
||||
s.p.is_finite()
|
||||
&& s.up.is_finite()
|
||||
&& s.across.is_none_or(|a| a.leave.is_finite() && a.normal.is_finite() && a.bend.is_none_or(|b| b.is_finite()))
|
||||
&& s.across.is_none_or(|a| {
|
||||
a.leave.is_finite() && a.normal.is_finite() && a.bend.is_none_or(|b| b.is_finite()) && a.jet.is_none_or(|j| j.is_finite())
|
||||
})
|
||||
}
|
||||
|
||||
/// How close a grid came: what [`fill_general`] reports when none is close enough.
|
||||
|
||||
@@ -30,7 +30,7 @@ fn tube_end(r: f64, g2: bool) -> impl FnMut(usize, &[f64]) -> Result<Vec<Station
|
||||
.iter()
|
||||
.map(|t| {
|
||||
let a = t * TAU;
|
||||
let across = Across { leave: Vec3::Z, normal: radial(a), bend: g2.then_some(Vec3::ZERO) };
|
||||
let across = Across { leave: Vec3::Z, normal: radial(a), bend: g2.then_some(Vec3::ZERO), jet: None };
|
||||
Station { p: radial(a) * r, up: around(a).cross(Vec3::Z), across: Some(across), corner: false, clash: false }
|
||||
})
|
||||
.collect())
|
||||
@@ -52,7 +52,7 @@ fn window(r: f64, a: f64, h: f64) -> impl FnMut(usize, &[f64]) -> Result<Vec<Sta
|
||||
};
|
||||
let p = radial(angle) * r + Vec3::Z * z;
|
||||
let up = tangent.cross(leave);
|
||||
Station { p, up, across: Some(Across { leave, normal: radial(angle), bend: None }), corner: t == 0.0, clash: false }
|
||||
Station { p, up, across: Some(Across { leave, normal: radial(angle), bend: None, jet: None }), corner: t == 0.0, clash: false }
|
||||
})
|
||||
.collect())
|
||||
}
|
||||
@@ -118,7 +118,13 @@ fn a_planar_polygon_is_its_plane_and_a_point_lifts_it() {
|
||||
.iter()
|
||||
.map(|t| {
|
||||
let p = Vec3::new(a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t, 0.0);
|
||||
Station { p, up: Vec3::Z, across: Some(Across { leave: inward, normal: Vec3::Z, bend: None }), corner: *t == 0.0, clash: false }
|
||||
Station {
|
||||
p,
|
||||
up: Vec3::Z,
|
||||
across: Some(Across { leave: inward, normal: Vec3::Z, bend: None, jet: None }),
|
||||
corner: *t == 0.0,
|
||||
clash: false,
|
||||
}
|
||||
})
|
||||
.collect())
|
||||
};
|
||||
|
||||
+462
-144
@@ -1,6 +1,6 @@
|
||||
//! Filling a closed loop of exact B-spline curves with tensor B-spline patches whose boundary
|
||||
//! rows *are* those curves, so that every edge lies on the fill exactly and its pcurve is a
|
||||
//! straight line in (u, v), affine in the edge's own parameter.
|
||||
//! Filling a closed loop of exact B-spline curves, rational or not, with tensor B-spline patches
|
||||
//! whose boundary rows *are* those curves, so that every edge lies on the fill exactly and its
|
||||
//! pcurve is a straight line in (u, v), affine in the edge's own parameter.
|
||||
//!
|
||||
//! The formulation:
|
||||
//!
|
||||
@@ -14,23 +14,42 @@
|
||||
//! edge `k` along `v = 0`, spoke `k` up `u = 1`, spoke `k − 1` back along `v = 1` and the
|
||||
//! second half of edge `k − 1` down `u = 0`. Each spoke is one row of the two patches it
|
||||
//! divides, the same control points in both.
|
||||
//! - **Rational edges** (a circle's arc, as every conic is held) are boundary rows of a rational
|
||||
//! patch: the fill works on homogeneous control points `(w·P, w)`. An edge's homogeneous points
|
||||
//! may be scaled by a constant without changing its curve or its parameter, so each edge is
|
||||
//! scaled for the weights at the loop's corners to agree; that is possible exactly when the
|
||||
//! product of the edges' end-to-start weight ratios round the loop is one (SMLib's conics have
|
||||
//! weight one at their ends), and refused otherwise. Inside, the weights are the Coons blend of
|
||||
//! the boundary's and stay fixed, so the surface `S = Σ R_ij P_ij`, `R_ij = w_ij N_i N_j / W`,
|
||||
//! is linear in the free points and everything below is a linear least-squares problem.
|
||||
//! - **One space.** Every side, each on `[0, 1]` by an affine map of its parameter (which changes
|
||||
//! neither the curve nor its pieces), is made compatible with every other exactly: split into
|
||||
//! Bézier pieces at every break of any, each piece raised to the highest degree (NURBS Book
|
||||
//! §5.5), joined again with each break repeated as often as that degree
|
||||
//! ([`crate::spline_space::embed`]). Every patch is a net in that space both ways.
|
||||
//! Bézier pieces at its breaks, each piece raised to the fill's degree (NURBS Book §5.5), joined
|
||||
//! with each break repeated as often as that degree, then carried into the common space
|
||||
//! ([`crate::spline_space::embed`]). The common space repeats every side's breaks as often as
|
||||
//! the degree; the breaks refinement adds are simple knots, so the fill is as smooth across them
|
||||
//! as its degree allows. The degree is the sides' highest, and at least one more than the
|
||||
//! highest continuity order n asked (G1 quadratic … G4 quintic), so the fill's own
|
||||
//! derivatives up to order n are continuous across the knots refinement adds.
|
||||
//! - **The start**: each patch is the bilinearly blended Coons patch of its sides (S. A. Coons,
|
||||
//! 1967), built on the Greville abscissae (a linear function is reproduced exactly there), the
|
||||
//! spokes starting as straight lines to the centre, the mean of the midpoints.
|
||||
//! - **The inside**: with the loop's rows held, the other control points (each patch's interior,
|
||||
//! the spokes' and the centre, each once) minimise, by least squares, the misfit to the
|
||||
//! conditions plus a vanishing multiple of the patches' thin-plate energy
|
||||
//! `∫∫ |S_uu|² + 2|S_uv|² + |S_vv|²` (which alone decides a G0 fill):
|
||||
//! - G1 at stations along an edge: `n·∂S/∂ν = 0` (`ν` across the edge into the fill, `n` the
|
||||
//! neighbour's normal), with weak rows `∂S/∂ν = s·w` (`w` the neighbour's leave direction,
|
||||
//! `s` the side's own speed) that keep the fill leaving into itself;
|
||||
//! - G2: `n·∂²S/∂ν² = (∂S/∂ν·w)²·κ`, `κ` the neighbour's normal curvature along `w`, the speed
|
||||
//! read from the last solution;
|
||||
//! `∫∫ |S_uu|² + 2|S_uv|² + |S_vv|²` (of the control net on the polynomial basis; it alone
|
||||
//! decides a G0 fill). The conditions are those of G^n continuity in the linkage-curve
|
||||
//! formulation ([`crate::linkage`]): at stations along an edge, the iso-parameter curve `c`
|
||||
//! of the fill leaving across it (σ growing into the fill) must have its height above the
|
||||
//! neighbour vanish to order n:
|
||||
//! - order 1: `n·c′ = 0` (`n` the neighbour's normal), with weak rows `w·c′ = s`, `t·c′ = 0`
|
||||
//! (`w` the neighbour's leave direction, `t = n × w` along the edge, `s` how fast the Coons
|
||||
//! start leaves along `w`) that keep the fill leaving into itself, about as fast; they act
|
||||
//! in the tangent plane only, so they never pull against the G1 row;
|
||||
//! - order k = 2 … n: `n·c⁽ᵏ⁾ = k!·[σᵏ] Σ_m H_m(x(σ))/m!` ([`crate::linkage::Jet::targets`]):
|
||||
//! the normal curvature (G2), its rate of change (G3) and that rate's (G4) are the
|
||||
//! neighbour's, in the fill's own leave direction. The right-hand side reads `c′ … c⁽ᵏ⁻¹⁾`
|
||||
//! from the last solution, so the rows are solved again, a few rounds;
|
||||
//! - across each spoke, the two patches' derivatives away from it opposed (`∂S_k/∂u + ∂S_{k+1}/∂v
|
||||
//! = 0` on the shared parameter): the patches join tangent there;
|
||||
//! - points the fill passes through, located on the last solution (nearest point of any patch)
|
||||
@@ -40,18 +59,20 @@
|
||||
//! the centre the spokes' tangency conditions overdetermine the patches: the stations near
|
||||
//! either are fitted but not held to the angle tolerance.
|
||||
//! - **Refinement**: while the normals miss the neighbours' (or each other's across a spoke) by
|
||||
//! more than the angle tolerance, a point is missed, or a patch folds, every span is halved
|
||||
//! (knot insertion carries the sides into the finer space exactly), up to [`MAX_ISO_SPANS`].
|
||||
//! more than the angle tolerance, a jump of order 2 … n carried over the loop's size
|
||||
//! ([`within_order`]) does, a point is missed, a patch folds, or the space has too few
|
||||
//! rows for the conditions from opposite sides to stay apart, every span is halved (knot
|
||||
//! insertion carries the sides into the finer space exactly), up to [`MAX_ISO_SPANS`].
|
||||
//! - **A crown** fills the middle out along the normals, on control points clear of the rows the
|
||||
//! conditions hold, so the loop and the conditions across it do not move.
|
||||
//!
|
||||
//! The fill's boundary is exact, whatever its inside: a station's gap is nought by construction.
|
||||
//! Rational edges are not supported (the Coons sum needs one weight function for all sides).
|
||||
//! Input is hostile: sizes are capped and every evaluation checked; nothing here panics.
|
||||
|
||||
use crate::Vec3;
|
||||
use crate::bspline::{HPoint, MAX_SPLINE_DEGREE, NurbsSurface, SplineError, SplineResult, check_knots, invalid, segment};
|
||||
use crate::fill_general::Across;
|
||||
use crate::linkage::MAX_CONTINUITY;
|
||||
use crate::spline_space::{Space, embed};
|
||||
|
||||
/// Most spans of the fill's space.
|
||||
@@ -62,17 +83,28 @@ pub const MAX_ISO_EDGES: usize = 16;
|
||||
const MAX_UNKNOWN_POINTS: usize = 500;
|
||||
/// Stations per span of a side where conditions are asked and checked.
|
||||
const STATIONS_PER_SPAN: usize = 6;
|
||||
/// Stations per control point of a row.
|
||||
const STATIONS_PER_POINT: usize = 3;
|
||||
/// Fewest stations along a side.
|
||||
const MIN_STATIONS: usize = 12;
|
||||
/// Weight of the weak Hermite rows against the G1 rows.
|
||||
/// Weight of the weak Hermite rows (in the tangent plane) against the G1 rows.
|
||||
const HERMITE_WEIGHT: f64 = 0.05;
|
||||
/// Rounds of relinearising G2 and relocating the points.
|
||||
/// Rounds of relinearising the rows of order two and up and relocating the points, per order
|
||||
/// above one (each order's right-hand side reads the lower derivatives of the last solution).
|
||||
const ROUNDS_PER_ORDER: usize = 2;
|
||||
/// Rounds when only points are relocated.
|
||||
const ROUNDS: usize = 3;
|
||||
/// The share of a side, from a clashing corner or the centre, whose stations are not held to
|
||||
/// the angle tolerance.
|
||||
const CLASH_SHARE: f64 = 0.15;
|
||||
/// The share of a side, from either of its corners, whose stations' mismatches of order two and
|
||||
/// up are not reported: at a corner the two sides' cross-boundary conditions meet and
|
||||
/// overdetermine the fill (the compatibility of boundary data at a vertex).
|
||||
const CORNER_SHARE: f64 = 0.1;
|
||||
/// How many points each way the fold check samples.
|
||||
const FOLD_SAMPLES: usize = 24;
|
||||
/// How far, relative to one, the corners' weights may disagree once the edges are scaled.
|
||||
const WEIGHT_TOL: f64 = 1e-9;
|
||||
/// Four-point Gauss–Legendre nodes and weights on `[−1, 1]` (exact to degree 7).
|
||||
const GAUSS4: [(f64, f64); 4] = [
|
||||
(-0.861_136_311_594_052_6, 0.347_854_845_137_453_8),
|
||||
@@ -81,13 +113,15 @@ const GAUSS4: [(f64, f64); 4] = [
|
||||
(0.861_136_311_594_052_6, 0.347_854_845_137_453_8),
|
||||
];
|
||||
|
||||
/// One edge of the loop: a non-rational B-spline from the loop's corner before it to the one
|
||||
/// after (the loop's direction).
|
||||
/// One edge of the loop: a B-spline from the loop's corner before it to the one after (the
|
||||
/// loop's direction), rational when it has weights.
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
pub struct IsoEdge {
|
||||
pub degree: usize,
|
||||
pub knots: Vec<f64>,
|
||||
pub points: Vec<Vec3>,
|
||||
/// One positive weight per point; `None` for a polynomial edge.
|
||||
pub weights: Option<Vec<f64>>,
|
||||
}
|
||||
|
||||
impl IsoEdge {
|
||||
@@ -101,6 +135,12 @@ impl IsoEdge {
|
||||
let [a, b] = self.domain();
|
||||
0.5 * (a + b)
|
||||
}
|
||||
|
||||
/// The homogeneous control points.
|
||||
fn homogeneous(&self) -> Vec<HPoint> {
|
||||
let w = |i: usize| self.weights.as_ref().and_then(|w| w.get(i)).copied().unwrap_or(1.0);
|
||||
self.points.iter().enumerate().map(|(i, p)| HPoint::weighted(*p, w(i))).collect()
|
||||
}
|
||||
}
|
||||
|
||||
/// What the fill is asked for besides its boundary.
|
||||
@@ -138,8 +178,8 @@ pub struct IsoPatch {
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
pub struct IsoFill {
|
||||
pub patches: Vec<IsoPatch>,
|
||||
/// Each spoke as a curve (degree, knots on `[0, 1]`, control points): from its edge's
|
||||
/// midpoint to the centre. Empty for a fill of one patch.
|
||||
/// Each spoke as a curve (degree, knots on `[0, 1]`, control points and, on a rational fill,
|
||||
/// weights): from its edge's midpoint to the centre. Empty for a fill of one patch.
|
||||
pub spokes: Vec<IsoEdge>,
|
||||
pub achieved: Achieved,
|
||||
}
|
||||
@@ -151,28 +191,30 @@ pub struct Achieved {
|
||||
/// a spoke, between the two patches') at a station held to the tolerance, and at any station.
|
||||
pub normal_deviation: f64,
|
||||
pub normal_deviation_all: f64,
|
||||
/// The largest difference (1/length) between the fill's normal curvature across a G2 edge
|
||||
/// and the neighbour's.
|
||||
pub curvature_mismatch: f64,
|
||||
/// `jumps[k]`, k = 2 … 4: the largest jump, clear of the loop's corners, of the k-th
|
||||
/// derivative of the linkage curves' height above the neighbour per unit length of their
|
||||
/// foot (`(n·c⁽ᵏ⁾ − target) / |x₁|ᵏ`, 1/lengthᵏ⁻¹) across the edges asked for order k or more:
|
||||
/// the normal curvature (k = 2), its rate of change (3) and that rate's (4). Zero below 2.
|
||||
pub jumps: [f64; MAX_CONTINUITY + 1],
|
||||
/// The farthest point from the fill.
|
||||
pub point_deviation: f64,
|
||||
}
|
||||
|
||||
/// A curve on `[0, 1]`.
|
||||
/// A curve on `[0, 1]`, on homogeneous control points.
|
||||
#[derive(Clone, Debug)]
|
||||
struct Curve {
|
||||
space: Space,
|
||||
points: Vec<Vec3>,
|
||||
points: Vec<HPoint>,
|
||||
}
|
||||
|
||||
/// The piece `[a, b]` of `edge`, on `[0, 1]` (run backwards when `reverse`).
|
||||
fn on_unit(edge: &IsoEdge, [a, b]: [f64; 2], reverse: bool) -> SplineResult<Curve> {
|
||||
/// The piece `[a, b]` of `edge` (homogeneous points scaled by `scale`), on `[0, 1]` (run
|
||||
/// backwards when `reverse`).
|
||||
fn on_unit(edge: &IsoEdge, scale: f64, [a, b]: [f64; 2], reverse: bool) -> SplineResult<Curve> {
|
||||
let p = edge.degree;
|
||||
let h: Vec<HPoint> = edge.points.iter().map(|q| HPoint::weighted(*q, 1.0)).collect();
|
||||
let (knots, h) = segment(p, &edge.knots, &h, a, b)?;
|
||||
let h: Vec<HPoint> = edge.homogeneous().into_iter().map(|q| q * scale).collect();
|
||||
let (knots, mut points) = segment(p, &edge.knots, &h, a, b)?;
|
||||
let map = |k: f64| if reverse { (b - k) / (b - a) } else { (k - a) / (b - a) };
|
||||
let mut knots: Vec<f64> = knots.iter().map(|k| map(*k)).collect();
|
||||
let mut points: Vec<Vec3> = h.iter().map(|q| q.xyz()).collect();
|
||||
if reverse {
|
||||
knots.reverse();
|
||||
points.reverse();
|
||||
@@ -188,8 +230,9 @@ fn on_unit(edge: &IsoEdge, [a, b]: [f64; 2], reverse: bool) -> SplineResult<Curv
|
||||
Ok(Curve { space: Space::new(p, knots)?, points })
|
||||
}
|
||||
|
||||
/// The straight line from `a` to `b` on `[0, 1]` (a spoke's start, or a pole when `a == b`).
|
||||
fn line(a: Vec3, b: Vec3) -> SplineResult<Curve> {
|
||||
/// The straight line from `a` to `b` on `[0, 1]`, homogeneous (a spoke's start, or a pole when
|
||||
/// the points are one).
|
||||
fn line(a: HPoint, b: HPoint) -> SplineResult<Curve> {
|
||||
Ok(Curve { space: Space::polynomial(1, 0.0, 1.0)?, points: vec![a, b] })
|
||||
}
|
||||
|
||||
@@ -200,32 +243,51 @@ fn breaks(s: &Space) -> Vec<f64> {
|
||||
out
|
||||
}
|
||||
|
||||
/// One space holding every curve: every break of any (and `extra`), each repeated as often as the
|
||||
/// highest degree; and each curve's control points there (each Bézier piece raised to that
|
||||
/// degree). Exact.
|
||||
fn compatible(curves: &[&Curve], extra: &[f64]) -> SplineResult<(Space, Vec<Vec<Vec3>>)> {
|
||||
let q = curves.iter().map(|s| s.space.degree()).max().unwrap_or(1);
|
||||
let mut all: Vec<f64> = curves.iter().flat_map(|s| breaks(&s.space)).chain(extra.iter().copied()).filter(|x| (0.0..=1.0).contains(x)).collect();
|
||||
all.sort_by(f64::total_cmp);
|
||||
all.dedup();
|
||||
/// Breaks of different sides closer than this (on `[0, 1]`) are one break: the rounding of
|
||||
/// mapping each side's parameter onto the unit interval.
|
||||
const BREAK_TOL: f64 = 1e-12;
|
||||
|
||||
/// The knots of `s` with each one within [`BREAK_TOL`] of one of `breaks` moved onto it.
|
||||
fn snapped(s: &Space, breaks: &[f64]) -> SplineResult<Space> {
|
||||
let knots = s.knots().iter().map(|k| breaks.iter().copied().find(|b| (b - k).abs() <= BREAK_TOL).unwrap_or(*k)).collect();
|
||||
Space::new(s.degree(), knots)
|
||||
}
|
||||
|
||||
/// One space holding every curve, of degree `q` (at least each curve's): every break of any
|
||||
/// curve repeated `q` times (breaks of different curves within [`BREAK_TOL`] made one), and the
|
||||
/// breaks `extra` adds as simple knots; and each curve's control points there (each Bézier piece
|
||||
/// raised to that degree, then knots inserted). Exact.
|
||||
fn compatible(curves: &[&Curve], extra: &[f64], q: usize) -> SplineResult<(Space, Vec<Vec<HPoint>>)> {
|
||||
let q = curves.iter().map(|s| s.space.degree()).max().unwrap_or(1).max(q);
|
||||
let mut own: Vec<f64> = curves.iter().flat_map(|s| breaks(&s.space)).collect();
|
||||
own.sort_by(f64::total_cmp);
|
||||
own.dedup_by(|b, a| (*b - *a).abs() <= BREAK_TOL);
|
||||
let inner: Vec<f64> = own.iter().copied().filter(|x| *x > BREAK_TOL && *x < 1.0 - BREAK_TOL).collect();
|
||||
let mut knots = vec![0.0; q + 1];
|
||||
for &x in all.iter().filter(|x| **x > 0.0 && **x < 1.0) {
|
||||
for &x in &inner {
|
||||
knots.extend(std::iter::repeat_n(x, q));
|
||||
}
|
||||
knots.extend(std::iter::repeat_n(1.0, q + 1));
|
||||
let target = Space::new(q, knots)?;
|
||||
let target = Space::new(q, knots)?.refined(extra)?;
|
||||
let mut out = Vec::with_capacity(curves.len());
|
||||
for c in curves {
|
||||
let p = c.space.degree();
|
||||
let h: Vec<HPoint> = c.points.iter().map(|x| HPoint::weighted(*x, 1.0)).collect();
|
||||
let mut points: Vec<Vec3> = Vec::new();
|
||||
for w in all.windows(2) {
|
||||
let space = snapped(&c.space, &[&[0.0], inner.as_slice(), &[1.0]].concat())?;
|
||||
let mine = breaks(&space);
|
||||
// The curve at degree q on its own breaks, each repeated q times.
|
||||
let mut points: Vec<HPoint> = Vec::new();
|
||||
for w in mine.windows(2) {
|
||||
let (a, b) = (w[0], w[1]);
|
||||
let (_, piece) = segment(p, c.space.knots(), &h, a, b)?;
|
||||
let piece: Vec<Vec3> = piece.iter().map(|x| x.xyz()).collect();
|
||||
let (_, piece) = segment(p, space.knots(), &c.points, a, b)?;
|
||||
let raised = embed(&Space::polynomial(p, a, b)?, &piece, &Space::polynomial(q, a, b)?)?;
|
||||
points.extend(raised.into_iter().skip(usize::from(!points.is_empty())));
|
||||
}
|
||||
let mut own_knots = vec![0.0; q + 1];
|
||||
for &x in mine.iter().filter(|x| **x > 0.0 && **x < 1.0) {
|
||||
own_knots.extend(std::iter::repeat_n(x, q));
|
||||
}
|
||||
own_knots.extend(std::iter::repeat_n(1.0, q + 1));
|
||||
let points = embed(&Space::new(q, own_knots)?, &points, &target)?;
|
||||
if points.len() != target.count() {
|
||||
return invalid("the sides' pieces do not fill the common space");
|
||||
}
|
||||
@@ -241,11 +303,14 @@ fn greville(s: &Space) -> Vec<f64> {
|
||||
}
|
||||
|
||||
/// A tensor patch in the fill's space both ways, control point `(i, j)` (i along u) at
|
||||
/// `i + c·j`, `c` the space's count.
|
||||
/// `i + c·j`, `c` the space's count, with its fixed weights.
|
||||
#[derive(Clone, Debug)]
|
||||
struct Patch {
|
||||
space: Space,
|
||||
net: Vec<Vec3>,
|
||||
weights: Vec<f64>,
|
||||
/// Whether any weight differs from one (else the polynomial basis is used as it is).
|
||||
rational: bool,
|
||||
/// The patch collapses along `u = 0` (a three-sided patch).
|
||||
pole: bool,
|
||||
}
|
||||
@@ -259,15 +324,86 @@ impl Patch {
|
||||
self.net.get(i + self.count() * j).copied().unwrap_or(Vec3::ZERO)
|
||||
}
|
||||
|
||||
fn weight(&self, i: usize, j: usize) -> f64 {
|
||||
self.weights.get(i + self.count() * j).copied().unwrap_or(1.0)
|
||||
}
|
||||
|
||||
/// The coefficients of `∂^(du+dv) S / ∂u^du ∂v^dv` at (u, v) on each control point it
|
||||
/// depends on: `(i, j, coefficient)`.
|
||||
/// depends on: `(i, j, coefficient)`. On a rational patch these are the derivatives of the
|
||||
/// rational basis `R_ij = w_ij N_i N_j / W`, by the Leibniz rule for `w_ij N_i N_j = R_ij W`:
|
||||
/// `R^(a,b) = (f^(a,b) − Σ_{(k,l) ≠ (0,0)} C(a,k) C(b,l) W^(k,l) R^(a−k,b−l)) / W`.
|
||||
fn row(&self, u: f64, v: f64, du: usize, dv: usize) -> SplineResult<Vec<(usize, usize, f64)>> {
|
||||
let (iu, bu) = self.space.basis(u, du)?;
|
||||
let (iv, bv) = self.space.basis(v, dv)?;
|
||||
let mut out = Vec::with_capacity(bu.len() * bv.len());
|
||||
for (b, nv) in bv.iter().enumerate() {
|
||||
for (a, nu) in bu.iter().enumerate() {
|
||||
out.push((iu + a, iv + b, nu * nv));
|
||||
if !self.rational {
|
||||
let (iu, bu) = self.space.basis(u, du)?;
|
||||
let (iv, bv) = self.space.basis(v, dv)?;
|
||||
let mut out = Vec::with_capacity(bu.len() * bv.len());
|
||||
for (b, nv) in bv.iter().enumerate() {
|
||||
for (a, nu) in bu.iter().enumerate() {
|
||||
out.push((iu + a, iv + b, nu * nv));
|
||||
}
|
||||
}
|
||||
return Ok(out);
|
||||
}
|
||||
// The basis and its derivatives up to the orders asked, each way.
|
||||
let ders = |t: f64, top: usize| -> SplineResult<(usize, Vec<Vec<f64>>)> {
|
||||
let mut first = 0;
|
||||
let mut out = Vec::with_capacity(top + 1);
|
||||
for r in 0..=top {
|
||||
let (f, b) = self.space.basis(t, r)?;
|
||||
first = f;
|
||||
out.push(b);
|
||||
}
|
||||
Ok((first, out))
|
||||
};
|
||||
let (iu, bu) = ders(u, du)?;
|
||||
let (iv, bv) = ders(v, dv)?;
|
||||
let width = (bu.first().map_or(0, Vec::len), bv.first().map_or(0, Vec::len));
|
||||
let basis = |d: &[Vec<f64>], r: usize, a: usize| d.get(r).and_then(|x| x.get(a)).copied().unwrap_or(0.0);
|
||||
// W^(k,l) at the point.
|
||||
let mut wd = vec![vec![0.0; dv + 1]; du + 1];
|
||||
for b in 0..width.1 {
|
||||
for a in 0..width.0 {
|
||||
let w = self.weight(iu + a, iv + b);
|
||||
for (k, row) in wd.iter_mut().enumerate() {
|
||||
for (l, x) in row.iter_mut().enumerate() {
|
||||
*x += w * basis(&bu, k, a) * basis(&bv, l, b);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
let w0 = wd.first().and_then(|r| r.first()).copied().unwrap_or(0.0);
|
||||
if !(w0 > 0.0) {
|
||||
return invalid("a patch's weights vanish");
|
||||
}
|
||||
let binom = |n: usize, k: usize| (0..k).fold(1.0, |c, i| c * (n - i) as f64 / (i + 1) as f64);
|
||||
let mut out = Vec::with_capacity(width.0 * width.1);
|
||||
let mut r = vec![vec![0.0; dv + 1]; du + 1];
|
||||
for b in 0..width.1 {
|
||||
for a in 0..width.0 {
|
||||
let w = self.weight(iu + a, iv + b);
|
||||
for total in 0..=du + dv {
|
||||
for k in 0..=du.min(total) {
|
||||
let l = total - k;
|
||||
if l > dv {
|
||||
continue;
|
||||
}
|
||||
let mut x = w * basis(&bu, k, a) * basis(&bv, l, b);
|
||||
for i in 0..=k {
|
||||
for j in 0..=l {
|
||||
if i == 0 && j == 0 {
|
||||
continue;
|
||||
}
|
||||
let wij = wd.get(i).and_then(|x| x.get(j)).copied().unwrap_or(0.0);
|
||||
let rr = r.get(k - i).and_then(|x| x.get(l - j)).copied().unwrap_or(0.0);
|
||||
x -= binom(k, i) * binom(l, j) * wij * rr;
|
||||
}
|
||||
}
|
||||
if let Some(slot) = r.get_mut(k).and_then(|x| x.get_mut(l)) {
|
||||
*slot = x / w0;
|
||||
}
|
||||
}
|
||||
}
|
||||
out.push((iu + a, iv + b, r.get(du).and_then(|x| x.get(dv)).copied().unwrap_or(0.0)));
|
||||
}
|
||||
}
|
||||
Ok(out)
|
||||
@@ -284,7 +420,8 @@ impl Patch {
|
||||
|
||||
fn surface(&self) -> SplineResult<NurbsSurface> {
|
||||
let c = self.count();
|
||||
let ctrl: Vec<HPoint> = (0..c).flat_map(|j| (0..c).map(move |i| (i, j))).map(|(i, j)| HPoint::weighted(self.at(i, j), 1.0)).collect();
|
||||
let ctrl: Vec<HPoint> =
|
||||
(0..c).flat_map(|j| (0..c).map(move |i| (i, j))).map(|(i, j)| HPoint::weighted(self.at(i, j), self.weight(i, j))).collect();
|
||||
let knots = self.space.knots().to_vec();
|
||||
NurbsSurface::from_homogeneous([self.space.degree(); 2], [knots.clone(), knots], ctrl)
|
||||
}
|
||||
@@ -327,19 +464,48 @@ struct Quilt {
|
||||
spokes: Vec<(usize, usize)>,
|
||||
}
|
||||
|
||||
/// The constants each edge's homogeneous points are scaled by so that the weights at every
|
||||
/// corner of the loop agree (module docs); refused when no such constants exist.
|
||||
fn corner_scales(edges: &[IsoEdge]) -> SplineResult<Vec<f64>> {
|
||||
let ends = |e: &IsoEdge| {
|
||||
let w = e.weights.as_ref();
|
||||
(w.and_then(|w| w.first()).copied().unwrap_or(1.0), w.and_then(|w| w.last()).copied().unwrap_or(1.0))
|
||||
};
|
||||
let mut scales = Vec::with_capacity(edges.len());
|
||||
let mut c = 1.0 / edges.first().map_or(1.0, |e| ends(e).0);
|
||||
for (k, e) in edges.iter().enumerate() {
|
||||
if k > 0 {
|
||||
let (_, before) = edges.get(k - 1).map(ends).unwrap_or((1.0, 1.0));
|
||||
c *= before / ends(e).0;
|
||||
}
|
||||
scales.push(c);
|
||||
}
|
||||
let closing = edges.last().map_or(1.0, |e| ends(e).1) * scales.last().copied().unwrap_or(1.0);
|
||||
if !(closing.is_finite() && (closing - 1.0).abs() <= WEIGHT_TOL) {
|
||||
return Err(SplineError::Unsupported(
|
||||
"a patch on rational edges whose weights at the loop's corners cannot be made to agree (their end-to-start ratios round the loop multiply to more or less than one)".into(),
|
||||
));
|
||||
}
|
||||
Ok(scales)
|
||||
}
|
||||
|
||||
impl Quilt {
|
||||
/// The fill's patches over the loop's sides in the space refined by `extra` breaks.
|
||||
fn new(edges: &[IsoEdge], extra: &[f64]) -> SplineResult<Quilt> {
|
||||
/// The fill's patches over the loop's sides at degree `degree` (at least) in the space
|
||||
/// refined by `extra` breaks.
|
||||
fn new(edges: &[IsoEdge], extra: &[f64], degree: usize) -> SplineResult<Quilt> {
|
||||
let n = edges.len();
|
||||
let corner = |k: usize| edges.get(k).and_then(|e| e.points.first()).copied().unwrap_or(Vec3::ZERO);
|
||||
let scales = corner_scales(edges)?;
|
||||
let scale = |k: usize| scales.get(k).copied().unwrap_or(1.0);
|
||||
let corner = |k: usize| edges.get(k).and_then(|e| e.homogeneous().first().copied()).map_or(HPoint::ZERO, |h| h * scale(k));
|
||||
if n <= 4 {
|
||||
// One patch: the edges are its sides (side 2 and 3 run back along their rows).
|
||||
let mut sides: Vec<Curve> = edges.iter().enumerate().map(|(k, e)| on_unit(e, e.domain(), k >= 2)).collect::<SplineResult<_>>()?;
|
||||
let mut sides: Vec<Curve> =
|
||||
edges.iter().enumerate().map(|(k, e)| on_unit(e, scale(k), e.domain(), k >= 2)).collect::<SplineResult<_>>()?;
|
||||
if n == 3 {
|
||||
sides.push(line(corner(0), corner(0))?);
|
||||
}
|
||||
let refs: Vec<&Curve> = sides.iter().collect();
|
||||
let (space, rows) = compatible(&refs, extra)?;
|
||||
let (space, rows) = compatible(&refs, extra, degree)?;
|
||||
let (Some(c0), Some(c1), Some(c2), Some(c3)) = (rows.first(), rows.get(1), rows.get(2), rows.get(3)) else { return invalid("sides") };
|
||||
let patch = coons(&space, c0, c2, c3, c1, n == 3)?;
|
||||
let c = space.count();
|
||||
@@ -357,20 +523,20 @@ impl Quilt {
|
||||
return Ok(Quilt { patches: vec![patch], slots: vec![slots], unknowns, pieces, spokes: Vec::new() });
|
||||
}
|
||||
// N patches about the centre: halves of each edge, and straight spokes to start.
|
||||
let middle: Vec<Vec3> = edges.iter().map(|e| eval_edge(e, e.middle())).collect::<SplineResult<_>>()?;
|
||||
let centre = middle.iter().fold(Vec3::ZERO, |s, p| s + *p) * (1.0 / n as f64);
|
||||
let middle: Vec<HPoint> = edges.iter().enumerate().map(|(k, e)| Ok(eval_edge_h(e, e.middle())? * scale(k))).collect::<SplineResult<_>>()?;
|
||||
let centre = middle.iter().filter_map(|h| h.project()).fold(Vec3::ZERO, |s, p| s + p) * (1.0 / n as f64);
|
||||
let mut curves: Vec<Curve> = Vec::with_capacity(3 * n);
|
||||
for e in edges {
|
||||
for (k, e) in edges.iter().enumerate() {
|
||||
let [a, b] = e.domain();
|
||||
curves.push(on_unit(e, [a, e.middle()], false)?);
|
||||
curves.push(on_unit(e, scale(k), [a, e.middle()], false)?);
|
||||
// The second half runs up u = 0 of the next patch: from V (v = 0) back to M.
|
||||
curves.push(on_unit(e, [e.middle(), b], true)?);
|
||||
curves.push(on_unit(e, scale(k), [e.middle(), b], true)?);
|
||||
}
|
||||
for m in &middle {
|
||||
curves.push(line(*m, centre)?);
|
||||
curves.push(line(*m, HPoint::weighted(centre, 1.0))?);
|
||||
}
|
||||
let refs: Vec<&Curve> = curves.iter().collect();
|
||||
let (space, rows) = compatible(&refs, extra)?;
|
||||
let (space, rows) = compatible(&refs, extra, degree)?;
|
||||
let c = space.count();
|
||||
let row = |r: usize| rows.get(r).cloned().ok_or_else(|| SplineError::Invalid("rows".into()));
|
||||
let mut patches = Vec::with_capacity(n);
|
||||
@@ -380,6 +546,8 @@ impl Quilt {
|
||||
// u = 1: spoke k; v = 1: spoke k − 1, from M_{k−1} to C.
|
||||
patches.push(coons(&space, &row(2 * k)?, &row(2 * n + prev)?, &row(2 * prev + 1)?, &row(2 * n + k)?, false)?);
|
||||
}
|
||||
// The spokes are shared rows: their weights the same in both patches (the Coons blends
|
||||
// agree on them), and so are their points once solved.
|
||||
// Unknowns: each spoke's points past its midpoint, the centre once, each patch's
|
||||
// interior.
|
||||
let mut unknowns = 0usize;
|
||||
@@ -430,60 +598,115 @@ impl Quilt {
|
||||
}
|
||||
}
|
||||
|
||||
/// A point of an edge.
|
||||
fn eval_edge(e: &IsoEdge, t: f64) -> SplineResult<Vec3> {
|
||||
Space::new(e.degree, e.knots.clone())?.eval(&e.points, t, 0)
|
||||
/// A homogeneous point of an edge.
|
||||
fn eval_edge_h(e: &IsoEdge, t: f64) -> SplineResult<HPoint> {
|
||||
Space::new(e.degree, e.knots.clone())?.eval(&e.homogeneous(), t, 0)
|
||||
}
|
||||
|
||||
/// The bilinearly blended Coons patch: `c0`, `c1` along u at v = 0 and 1, `d0`, `d1` along v at
|
||||
/// u = 0 and 1, all running with their parameter.
|
||||
fn coons(space: &Space, c0: &[Vec3], c1: &[Vec3], d0: &[Vec3], d1: &[Vec3], pole: bool) -> SplineResult<Patch> {
|
||||
/// The bilinearly blended Coons patch of homogeneous sides: `c0`, `c1` along u at v = 0 and 1,
|
||||
/// `d0`, `d1` along v at u = 0 and 1, all running with their parameter. Points and weights are
|
||||
/// blended alike, so the boundary rows are the sides' own and the weights inside are the blend
|
||||
/// of theirs (refused where one is not positive).
|
||||
fn coons(space: &Space, c0: &[HPoint], c1: &[HPoint], d0: &[HPoint], d1: &[HPoint], pole: bool) -> SplineResult<Patch> {
|
||||
let g = greville(space);
|
||||
let c = space.count();
|
||||
let (Some(p00), Some(p10), Some(p01), Some(p11)) = (c0.first().copied(), c0.last().copied(), c1.first().copied(), c1.last().copied()) else {
|
||||
return invalid("a side has no points");
|
||||
};
|
||||
let get = |a: &[Vec3], k: usize| a.get(k).copied().unwrap_or(Vec3::ZERO);
|
||||
let get = |a: &[HPoint], k: usize| a.get(k).copied().unwrap_or(HPoint::ZERO);
|
||||
let mut net = Vec::with_capacity(c * c);
|
||||
let mut weights = Vec::with_capacity(c * c);
|
||||
for (j, &y) in g.iter().enumerate() {
|
||||
for (i, &x) in g.iter().enumerate() {
|
||||
let ruled_v = get(c0, i) * (1.0 - y) + get(c1, i) * y;
|
||||
let ruled_u = get(d0, j) * (1.0 - x) + get(d1, j) * x;
|
||||
let bilinear = (p00 * (1.0 - x) + p10 * x) * (1.0 - y) + (p01 * (1.0 - x) + p11 * x) * y;
|
||||
// The boundary rows exactly the sides' own points.
|
||||
net.push(match (i, j) {
|
||||
let h = match (i, j) {
|
||||
(_, 0) => get(c0, i),
|
||||
(_, j) if j + 1 == c => get(c1, i),
|
||||
(0, _) => get(d0, j),
|
||||
(i, _) if i + 1 == c => get(d1, j),
|
||||
_ => ruled_v + ruled_u - bilinear,
|
||||
});
|
||||
};
|
||||
let p = h
|
||||
.project()
|
||||
.filter(|_| h.w > 0.0)
|
||||
.ok_or_else(|| SplineError::Unsupported("a rational patch whose blended weights are not positive".into()))?;
|
||||
net.push(p);
|
||||
weights.push(h.w);
|
||||
}
|
||||
}
|
||||
Ok(Patch { space: space.clone(), net, pole })
|
||||
let rational = weights.iter().any(|w| *w != 1.0);
|
||||
Ok(Patch { space: space.clone(), net, weights, rational, pole })
|
||||
}
|
||||
|
||||
/// One condition row: coefficients on patches' control points `(patch, i, j, coefficient)`,
|
||||
/// dotted with `dir` when given, else on one coordinate `coord`; its value and weight.
|
||||
/// One condition row: `Σ along·P = value` over its terms' control points, and its weight.
|
||||
struct Row {
|
||||
terms: Vec<(usize, usize, usize, f64)>,
|
||||
dir: Option<Vec3>,
|
||||
coord: usize,
|
||||
terms: Vec<Term>,
|
||||
value: f64,
|
||||
weight: f64,
|
||||
}
|
||||
|
||||
/// A control point `(i, j)` of patch `patch` in a row, its coefficient a vector (dotted with
|
||||
/// the point).
|
||||
#[derive(Clone, Copy)]
|
||||
struct Term {
|
||||
patch: usize,
|
||||
i: usize,
|
||||
j: usize,
|
||||
along: Vec3,
|
||||
}
|
||||
|
||||
/// The terms of a patch's derivative coefficients `(i, j, c)` taken along `dir`.
|
||||
fn along(patch: usize, coefficients: &[(usize, usize, f64)], dir: Vec3) -> Vec<Term> {
|
||||
coefficients.iter().map(|&(i, j, c)| Term { patch, i, j, along: dir * c }).collect()
|
||||
}
|
||||
|
||||
/// The three rows `Σ terms·P = value` on each coordinate.
|
||||
fn on_each_axis(rows: &mut Vec<Row>, patch_terms: &[(usize, Vec<(usize, usize, f64)>)], value: Vec3, weight: f64) {
|
||||
for (axis, v) in [(Vec3::X, value.x), (Vec3::Y, value.y), (Vec3::Z, value.z)] {
|
||||
let terms = patch_terms.iter().flat_map(|(p, c)| along(*p, c, axis)).collect();
|
||||
rows.push(Row { terms, value: v, weight });
|
||||
}
|
||||
}
|
||||
|
||||
/// A station of a loop edge: its patch, side and fraction there, its condition, whether it is
|
||||
/// held to the angle tolerance, and the side's speed there.
|
||||
/// held to the angle tolerance, whether it is clear of the loop's corners (its curvature and
|
||||
/// curvature change are reported), and how fast the fill should leave there.
|
||||
struct Station {
|
||||
patch: usize,
|
||||
side: usize,
|
||||
t: f64,
|
||||
across: Across,
|
||||
held: bool,
|
||||
clear: bool,
|
||||
speed: f64,
|
||||
}
|
||||
|
||||
impl Station {
|
||||
/// The continuity order its conditions ask: 1 (G1), or its neighbour's jet's (G2 … G4).
|
||||
fn order(&self) -> usize {
|
||||
self.across.jet.map_or(1, |j| j.order)
|
||||
}
|
||||
}
|
||||
|
||||
/// The derivatives `[c′, c″, …, c⁽ⁿ⁾]` (n = `order`) of the linkage curve `c` leaving a station
|
||||
/// across its edge into the patch: the iso-parameter curve, σ growing into the patch.
|
||||
fn linkage(patch: &Patch, s: &Station, order: usize) -> SplineResult<Vec<Vec3>> {
|
||||
let (uv, (du, dv, sign)) = station_on(s.side, s.t);
|
||||
(1..=order).map(|k| Ok(patch.eval(uv[0], uv[1], k * du, k * dv)? * sign.powi(k as i32))).collect()
|
||||
}
|
||||
|
||||
/// What orders 2 … n ask of `n·c⁽ᵏ⁾` at a station ([`crate::linkage::Jet::targets`]), read
|
||||
/// from the linkage curve's last derivatives `c` (element k − 2 is order k's); empty for G1.
|
||||
fn targets(s: &Station, c: &[Vec3]) -> Vec<f64> {
|
||||
let Some(jet) = s.across.jet else { return Vec::new() };
|
||||
let (t, w) = (s.across.normal.cross(s.across.leave), s.across.leave);
|
||||
let x: Vec<[f64; 2]> = c.iter().map(|d| [d.dot(t), d.dot(w)]).collect();
|
||||
jet.targets(&x).unwrap_or_default()
|
||||
}
|
||||
|
||||
/// Fill the loop `edges` (three to [`MAX_ISO_EDGES`], each ending where the next starts) with
|
||||
/// patches on them (module docs). `across(k, ts)` gives the neighbour's conditions at the
|
||||
/// parameters `ts` of edge `k` (its own parameter; `None`: G0 there); `clash[k]`: the
|
||||
@@ -517,6 +740,11 @@ pub fn iso_fill(
|
||||
if !(a < b) || e.points.iter().any(|q| !q.is_finite()) {
|
||||
return invalid(format!("edge {k} has an empty domain or a control point that is not finite"));
|
||||
}
|
||||
if let Some(w) = &e.weights
|
||||
&& (w.len() != e.points.len() || w.iter().any(|x| !(x.is_finite() && *x > 0.0)))
|
||||
{
|
||||
return invalid(format!("edge {k}'s weights must be positive, one per control point"));
|
||||
}
|
||||
}
|
||||
let mut all = edges.iter().flat_map(|e| e.points.iter().copied());
|
||||
let first = all.next().ok_or_else(|| SplineError::Invalid("an edge has no control points".into()))?;
|
||||
@@ -529,19 +757,47 @@ pub fn iso_fill(
|
||||
return invalid(format!("edge {k} does not end where edge {} starts", (k + 1) % n));
|
||||
}
|
||||
}
|
||||
// The highest continuity asked anywhere, from the conditions at each edge's middle (an
|
||||
// edge's continuity is one along it).
|
||||
let mut order = 0usize;
|
||||
for k in 0..n {
|
||||
let Some(e) = edges.get(k) else { continue };
|
||||
let found = across(k, &[e.middle()]).map_err(SplineError::Invalid)?;
|
||||
for c in found.into_iter().flatten() {
|
||||
order = order.max(c.jet.map_or(1, |j| j.order));
|
||||
}
|
||||
}
|
||||
if order > MAX_CONTINUITY {
|
||||
return invalid(format!("a patch's continuity is at most G{MAX_CONTINUITY}"));
|
||||
}
|
||||
let degree = if order == 0 { 1 } else { (order + 1).min(MAX_SPLINE_DEGREE) };
|
||||
let mut extra: Vec<f64> = Vec::new();
|
||||
let mut spans = 1usize;
|
||||
loop {
|
||||
let quilt = Quilt::new(edges, &extra)?;
|
||||
let quilt = Quilt::new(edges, &extra, degree)?;
|
||||
let stations = stations(edges, &quilt, &mut *across, clash)?;
|
||||
let (quilt, report) = solve(quilt, &stations, spec)?;
|
||||
let fits = report.normal_deviation <= spec.angle_tol && report.point_deviation <= spec.tol && !quilt.patches.iter().any(folds);
|
||||
let space = quilt.patches.first().map(|p| p.space.clone()).ok_or_else(|| SplineError::Invalid("no patch".into()))?;
|
||||
// Rows enough for the conditions from opposite sides not to share control points.
|
||||
let roomy = space.count() >= 2 * (order + 1);
|
||||
let finer = quilt.unknowns * 4 > MAX_UNKNOWN_POINTS;
|
||||
if !roomy && !finer && spans * 2 <= MAX_ISO_SPANS {
|
||||
spans *= 2;
|
||||
let b = breaks(&space);
|
||||
extra.extend(b.windows(2).map(|w| 0.5 * (w[0] + w[1])));
|
||||
continue;
|
||||
}
|
||||
let (quilt, report) = solve(quilt, &stations, spec)?;
|
||||
let fits = report.normal_deviation <= spec.angle_tol
|
||||
&& within_order(&report.jumps, order, size, spec.angle_tol)
|
||||
&& report.point_deviation <= spec.tol
|
||||
&& !quilt.patches.iter().any(folds);
|
||||
if fits || spans * 2 > MAX_ISO_SPANS || finer {
|
||||
if !fits {
|
||||
let jumps =
|
||||
report.jumps.get(2..=order).filter(|_| order >= 2).map(|j| format!(", jumps of orders 2 to {order} {j:?}")).unwrap_or_default();
|
||||
let folded = if quilt.patches.iter().any(folds) { ", and it folds" } else { "" };
|
||||
return Err(SplineError::Unsupported(format!(
|
||||
"a patch on these edges within the tolerances: normals within {:.4}° of the neighbours', points within {:.2e}",
|
||||
"a patch on these edges within the tolerances: normals within {:.4}° of the neighbours'{jumps}, points within {:.2e}{folded}",
|
||||
report.normal_deviation.to_degrees(),
|
||||
report.point_deviation
|
||||
)));
|
||||
@@ -556,6 +812,21 @@ pub fn iso_fill(
|
||||
}
|
||||
}
|
||||
|
||||
/// Whether every jump of order 2 … `order` is within the angle tolerance once carried into the
|
||||
/// fill: by Taylor's theorem a jump `j` of order k (1/lengthᵏ⁻¹) turns the normal by
|
||||
/// `j·dᵏ⁻¹/(k − 1)!` at a distance `d` from the edge, and no point of the fill is farther from
|
||||
/// the loop than half its `size`. So each order is held as tightly as the tangent plane is.
|
||||
pub fn within_order(jumps: &[f64; MAX_CONTINUITY + 1], order: usize, size: f64, angle_tol: f64) -> bool {
|
||||
(2..=order).all(|k| jumps.get(k).is_some_and(|j| turn(*j, k, size) <= angle_tol))
|
||||
}
|
||||
|
||||
/// How far (radians) a jump `jump` of order `k` turns the normal over half of `size`
|
||||
/// ([`within_order`]).
|
||||
pub fn turn(jump: f64, k: usize, size: f64) -> f64 {
|
||||
let factorial: f64 = (1..k).map(|i| i as f64).product();
|
||||
jump * (0.5 * size).powi(k as i32 - 1) / factorial
|
||||
}
|
||||
|
||||
/// The patches, their sides and the spokes.
|
||||
fn finished(quilt: Quilt, edges: &[IsoEdge], achieved: Achieved) -> SplineResult<IsoFill> {
|
||||
let n = edges.len();
|
||||
@@ -585,12 +856,20 @@ fn finished(quilt: Quilt, edges: &[IsoEdge], achieved: Achieved) -> SplineResult
|
||||
.filter_map(|&(left, _)| {
|
||||
let p = quilt.patches.get(left)?;
|
||||
let c = p.count();
|
||||
Some(IsoEdge { degree: p.space.degree(), knots: p.space.knots().to_vec(), points: (0..c).map(|j| p.at(c - 1, j)).collect() })
|
||||
let weights = p.rational.then(|| (0..c).map(|j| p.weight(c - 1, j)).collect());
|
||||
Some(IsoEdge { degree: p.space.degree(), knots: p.space.knots().to_vec(), points: (0..c).map(|j| p.at(c - 1, j)).collect(), weights })
|
||||
})
|
||||
.collect();
|
||||
Ok(IsoFill { patches, spokes, achieved })
|
||||
}
|
||||
|
||||
/// Stations along a side: [`STATIONS_PER_SPAN`] a span and [`STATIONS_PER_POINT`] per control
|
||||
/// point of a row (each order of continuity holds one row of them, so the stations must outnumber
|
||||
/// a row's points well for the conditions to hold between them too), at least [`MIN_STATIONS`].
|
||||
fn station_count(space: &Space) -> usize {
|
||||
(breaks(space).len().saturating_sub(1) * STATIONS_PER_SPAN).max(STATIONS_PER_POINT * space.count()).max(MIN_STATIONS)
|
||||
}
|
||||
|
||||
/// The stations of every loop edge with a condition: enough of them for the space's spans.
|
||||
fn stations(
|
||||
edges: &[IsoEdge],
|
||||
@@ -601,7 +880,7 @@ fn stations(
|
||||
let n = edges.len();
|
||||
let mut out = Vec::new();
|
||||
let Some(space) = quilt.patches.first().map(|p| &p.space) else { return Ok(out) };
|
||||
let count = (breaks(space).len().saturating_sub(1) * STATIONS_PER_SPAN).max(MIN_STATIONS);
|
||||
let count = station_count(space);
|
||||
for &(k, [a, b], patch, side) in &quilt.pieces {
|
||||
let [ea, eb] = edges.get(k).map(IsoEdge::domain).unwrap_or([0.0, 1.0]);
|
||||
// Stations inside the piece (its corners are its neighbours' corners); every side runs
|
||||
@@ -616,9 +895,14 @@ fn stations(
|
||||
let f = (param - ea) / (eb - ea);
|
||||
let near_start = clash.get(k).copied().unwrap_or(false) && f < CLASH_SHARE;
|
||||
let near_end = clash.get((k + 1) % n).copied().unwrap_or(false) && f > 1.0 - CLASH_SHARE;
|
||||
let (uv, (du, dv, _)) = station_on(side, t);
|
||||
let speed = p.eval(uv[0], uv[1], dv, du).map(|d| d.len()).unwrap_or(0.0);
|
||||
out.push(Station { patch, side, t, across: c, held: !(near_start || near_end), speed });
|
||||
let clear = f > CORNER_SHARE && f < 1.0 - CORNER_SHARE;
|
||||
let (uv, (du, dv, sign)) = station_on(side, t);
|
||||
// How fast the Coons start leaves into the leave direction: it blends the sides
|
||||
// across, so it leaves as fast as the opposite sides lie apart; the side's own speed
|
||||
// where it does not leave into the patch at all.
|
||||
let along = p.eval(uv[0], uv[1], dv, du).map(|d| d.len()).unwrap_or(0.0);
|
||||
let speed = p.eval(uv[0], uv[1], du, dv).map(|d| (d * sign).dot(c.leave)).ok().filter(|x| *x > 0.0).unwrap_or(along);
|
||||
out.push(Station { patch, side, t, across: c, held: !(near_start || near_end), clear, speed });
|
||||
}
|
||||
}
|
||||
Ok(out)
|
||||
@@ -635,49 +919,34 @@ fn solve(mut quilt: Quilt, stations: &[Station], spec: &IsoSpec) -> SplineResult
|
||||
let energy: Vec<Vec<(usize, usize, f64)>> = quilt.patches.iter().map(energy).collect::<SplineResult<_>>()?;
|
||||
let mut located: Vec<(usize, [f64; 2])> = Vec::new();
|
||||
let mut report = Achieved::default();
|
||||
let spoke_count = quilt.patches.first().map_or(0, |p| (breaks(&p.space).len().saturating_sub(1) * STATIONS_PER_SPAN).max(MIN_STATIONS));
|
||||
for round in 0..ROUNDS {
|
||||
let spoke_count = quilt.patches.first().map_or(0, |p| station_count(&p.space));
|
||||
let top = stations.iter().map(Station::order).max().unwrap_or(1);
|
||||
// The space's span (the same both ways and in every patch).
|
||||
let span = quilt.patches.first().map_or(1.0, |p| 1.0 / breaks(&p.space).len().saturating_sub(1).max(1) as f64);
|
||||
let rounds = if top >= 2 { ROUNDS_PER_ORDER * top } else { ROUNDS };
|
||||
for round in 0..rounds {
|
||||
let mut rows: Vec<Row> = Vec::new();
|
||||
for s in stations {
|
||||
let Some(patch) = quilt.patches.get(s.patch) else { continue };
|
||||
let (uv, (du, dv, sign)) = station_on(s.side, s.t);
|
||||
let first: Vec<(usize, usize, usize, f64)> =
|
||||
patch.row(uv[0], uv[1], du, dv)?.into_iter().map(|(i, j, c)| (s.patch, i, j, c * sign)).collect();
|
||||
rows.push(Row { terms: first.clone(), dir: Some(s.across.normal), coord: 0, value: 0.0, weight: 1.0 });
|
||||
let w = s.across.leave * s.speed;
|
||||
for (c, value) in [w.x, w.y, w.z].into_iter().enumerate() {
|
||||
rows.push(Row { terms: first.clone(), dir: None, coord: c, value, weight: HERMITE_WEIGHT });
|
||||
}
|
||||
if let Some(bend) = s.across.bend {
|
||||
let along = (patch.eval(uv[0], uv[1], du, dv)? * sign).dot(s.across.leave);
|
||||
let second = patch.row(uv[0], uv[1], 2 * du, 2 * dv)?.into_iter().map(|(i, j, c)| (s.patch, i, j, c)).collect();
|
||||
rows.push(Row { terms: second, dir: Some(s.across.normal), coord: 0, value: along * along * bend.dot(s.across.normal), weight: 1.0 });
|
||||
}
|
||||
linkage_rows(&mut rows, patch, s, span)?;
|
||||
}
|
||||
// Across each spoke the two patches' derivatives away from it are opposed.
|
||||
for &(left, right) in &quilt.spokes {
|
||||
let (Some(a), Some(b)) = (quilt.patches.get(left), quilt.patches.get(right)) else { continue };
|
||||
for i in 0..spoke_count {
|
||||
let t = (i as f64 + 0.5) / spoke_count as f64;
|
||||
let mut terms: Vec<(usize, usize, usize, f64)> = a.row(1.0, t, 1, 0)?.into_iter().map(|(i, j, c)| (left, i, j, c)).collect();
|
||||
terms.extend(b.row(t, 1.0, 0, 1)?.into_iter().map(|(i, j, c)| (right, i, j, c)));
|
||||
for c in 0..3 {
|
||||
rows.push(Row { terms: terms.clone(), dir: None, coord: c, value: 0.0, weight: 1.0 });
|
||||
}
|
||||
on_each_axis(&mut rows, &[(left, a.row(1.0, t, 1, 0)?), (right, b.row(t, 1.0, 0, 1)?)], Vec3::ZERO, 1.0);
|
||||
}
|
||||
}
|
||||
for (p, (patch, uv)) in spec.points.iter().zip(&located) {
|
||||
let Some(on) = quilt.patches.get(*patch) else { continue };
|
||||
let terms: Vec<(usize, usize, usize, f64)> = on.row(uv[0], uv[1], 0, 0)?.into_iter().map(|(i, j, c)| (*patch, i, j, c)).collect();
|
||||
for (c, value) in [p.x, p.y, p.z].into_iter().enumerate() {
|
||||
rows.push(Row { terms: terms.clone(), dir: None, coord: c, value, weight: 1.0 });
|
||||
}
|
||||
on_each_axis(&mut rows, &[(*patch, on.row(uv[0], uv[1], 0, 0)?)], *p, 1.0);
|
||||
}
|
||||
let solved = least_squares(&quilt, &rows, &energy)?;
|
||||
quilt.set(&solved);
|
||||
located = spec.points.iter().map(|p| nearest(&quilt, *p)).collect();
|
||||
report = measure(&quilt, stations, spec.points, &located, spoke_count)?;
|
||||
let relinearise = stations.iter().any(|s| s.across.bend.is_some()) || !spec.points.is_empty();
|
||||
let relinearise = top >= 2 || !spec.points.is_empty();
|
||||
if !relinearise && round == 0 {
|
||||
break;
|
||||
}
|
||||
@@ -685,6 +954,50 @@ fn solve(mut quilt: Quilt, stations: &[Station], spec: &IsoSpec) -> SplineResult
|
||||
Ok((quilt, report))
|
||||
}
|
||||
|
||||
/// A station's rows (module docs), linearised about the patch as it stands:
|
||||
/// - G1: `n·c′ = 0`, and the weak Hermite rows `w·c′ = s`, `t·c′ = 0` in the neighbour's tangent
|
||||
/// plane only (`t = n × w`), so they never pull against the G1 row;
|
||||
/// - orders k = 2 … n: `n·c⁽ᵏ⁾ = T_k(x₁ … x_{k−1})` ([`crate::linkage::Jet::linearised`]) by
|
||||
/// Newton's method: `T_k` is replaced by its tangent at the last solution's `x⁰`, so the row is
|
||||
/// `n·c⁽ᵏ⁾ − Σ_j ∇_j T_k · x_j = T_k(x⁰) − Σ_j ∇_j T_k · x_j⁰`, with `x_j = (t·c⁽ʲ⁾, w·c⁽ʲ⁾)`
|
||||
/// and `c⁽ʲ⁾ = signʲ ∂ʲS` across. Each is weighed by `hᵏ⁻¹` against the G1 row (`h` the span),
|
||||
/// so every order says how far it misses over one span and none swamps the others as the
|
||||
/// spans shrink.
|
||||
fn linkage_rows(rows: &mut Vec<Row>, patch: &Patch, s: &Station, span: f64) -> SplineResult<()> {
|
||||
let (uv, (du, dv, sign)) = station_on(s.side, s.t);
|
||||
let (n, w) = (s.across.normal, s.across.leave);
|
||||
let t = n.cross(w);
|
||||
let order = s.order();
|
||||
// The coefficients of c′ … c⁽ⁿ⁾ on the control points.
|
||||
let derivative: Vec<Vec<(usize, usize, f64)>> = (1..=order)
|
||||
.map(|k| {
|
||||
let sign_k = sign.powi(k as i32);
|
||||
Ok(patch.row(uv[0], uv[1], k * du, k * dv)?.into_iter().map(|(i, j, c)| (i, j, c * sign_k)).collect())
|
||||
})
|
||||
.collect::<SplineResult<_>>()?;
|
||||
let Some(first) = derivative.first() else { return Ok(()) };
|
||||
rows.push(Row { terms: along(s.patch, first, n), value: 0.0, weight: 1.0 });
|
||||
for (dir, value) in [(w, s.speed), (t, 0.0)] {
|
||||
rows.push(Row { terms: along(s.patch, first, dir), value, weight: HERMITE_WEIGHT });
|
||||
}
|
||||
let Some(jet) = s.across.jet else { return Ok(()) };
|
||||
let c = linkage(patch, s, order)?;
|
||||
let x0: Vec<[f64; 2]> = c.iter().map(|d| [d.dot(t), d.dot(w)]).collect();
|
||||
let Some(targets) = jet.linearised(&x0) else { return Ok(()) };
|
||||
for (k, target) in (2..).zip(targets) {
|
||||
let Some(top) = derivative.get(k - 1) else { continue };
|
||||
let mut terms = along(s.patch, top, n);
|
||||
let mut value = target.value;
|
||||
for (j, (g, x)) in target.gradient.iter().zip(&x0).enumerate().take(k - 1) {
|
||||
let Some(cj) = derivative.get(j) else { continue };
|
||||
terms.extend(along(s.patch, cj, (t * g[0] + w * g[1]) * -1.0));
|
||||
value -= g[0] * x[0] + g[1] * x[1];
|
||||
}
|
||||
rows.push(Row { terms, value, weight: span.powi(k as i32 - 1) });
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// A patch's thin-plate energy coefficients between its control points (flat indices), of the
|
||||
/// terms `G⁽ʳ⁾ ⊗ G⁽ˢ⁾` with `(r, s)` = (2, 0), (1, 1) twice, (0, 2).
|
||||
fn energy(patch: &Patch) -> SplineResult<Vec<(usize, usize, f64)>> {
|
||||
@@ -738,20 +1051,16 @@ fn least_squares(quilt: &Quilt, rows: &[Row], energy: &[Vec<(usize, usize, f64)>
|
||||
for r in rows {
|
||||
let mut unknown: Vec<(usize, f64)> = Vec::new();
|
||||
let mut known = 0.0;
|
||||
for &(p, i, j, c) in &r.terms {
|
||||
let flat = i + quilt.patches.get(p).map_or(0, Patch::count) * j;
|
||||
let coords = match r.dir {
|
||||
Some(d) => [d.x, d.y, d.z],
|
||||
None => std::array::from_fn(|k| if k == r.coord { 1.0 } else { 0.0 }),
|
||||
};
|
||||
let v = value(p, flat);
|
||||
for (k, f) in coords.iter().enumerate() {
|
||||
if *f == 0.0 {
|
||||
for t in &r.terms {
|
||||
let flat = t.i + quilt.patches.get(t.patch).map_or(0, Patch::count) * t.j;
|
||||
let v = value(t.patch, flat);
|
||||
for (k, (f, x)) in [t.along.x, t.along.y, t.along.z].into_iter().zip([v.x, v.y, v.z]).enumerate() {
|
||||
if f == 0.0 {
|
||||
continue;
|
||||
}
|
||||
match slot(p, flat) {
|
||||
Some(x) => unknown.push((3 * x + k, c * f)),
|
||||
None => known += c * f * [v.x, v.y, v.z][k],
|
||||
match slot(t.patch, flat) {
|
||||
Some(id) => unknown.push((3 * id + k, f)),
|
||||
None => known += f * x,
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -859,17 +1168,23 @@ fn measure(quilt: &Quilt, stations: &[Station], points: &[Vec3], located: &[(usi
|
||||
r.normal_deviation = r.normal_deviation.max(a);
|
||||
}
|
||||
};
|
||||
let mut mismatch = 0.0f64;
|
||||
let mut jumps = [0.0f64; MAX_CONTINUITY + 1];
|
||||
for s in stations {
|
||||
let Some(patch) = quilt.patches.get(s.patch) else { continue };
|
||||
let (uv, (du, dv, sign)) = station_on(s.side, s.t);
|
||||
let (uv, _) = station_on(s.side, s.t);
|
||||
let Some(normal) = patch.normal(uv[0], uv[1]) else { continue };
|
||||
angle(normal.dot(s.across.normal).abs().clamp(0.0, 1.0).acos(), s.held);
|
||||
if let Some(bend) = s.across.bend {
|
||||
let first = patch.eval(uv[0], uv[1], du, dv)? * sign;
|
||||
let second = patch.eval(uv[0], uv[1], 2 * du, 2 * dv)?;
|
||||
let mine = second.dot(s.across.normal) / first.len2().max(1e-300);
|
||||
mismatch = mismatch.max((mine - bend.dot(s.across.normal)).abs());
|
||||
if s.order() < 2 || !s.clear {
|
||||
continue;
|
||||
}
|
||||
// The linkage curve's height above the neighbour, per unit length along its foot.
|
||||
let c = linkage(patch, s, s.order())?;
|
||||
let n = s.across.normal;
|
||||
let Some(first) = c.first() else { continue };
|
||||
let foot = (*first - n * first.dot(n)).len().max(1e-300);
|
||||
for (k, want) in (2..).zip(targets(s, &c)) {
|
||||
let (Some(ck), Some(jump)) = (c.get(k - 1), jumps.get_mut(k)) else { continue };
|
||||
*jump = jump.max(((n.dot(*ck) - want) / foot.powi(k as i32)).abs());
|
||||
}
|
||||
}
|
||||
// Across the spokes, clear of the centre.
|
||||
@@ -882,7 +1197,7 @@ fn measure(quilt: &Quilt, stations: &[Station], points: &[Vec3], located: &[(usi
|
||||
}
|
||||
}
|
||||
}
|
||||
r.curvature_mismatch = mismatch;
|
||||
r.jumps = jumps;
|
||||
for (p, (patch, uv)) in points.iter().zip(located) {
|
||||
if let Some(on) = quilt.patches.get(*patch) {
|
||||
r.point_deviation = r.point_deviation.max(on.eval(uv[0], uv[1], 0, 0)?.dist(*p));
|
||||
@@ -946,15 +1261,14 @@ fn folds(patch: &Patch) -> bool {
|
||||
}
|
||||
|
||||
/// The fill with its middle filled out by `crown` times the loop's `size` along the normals, on
|
||||
/// control points clear of the rows the conditions hold (one row from a G0 edge, two from G1,
|
||||
/// three from G2), so the loop and the conditions across it do not move. One patch rises by
|
||||
/// `16·u(1−u)·v(1−v)`; patches about a centre by `u(2−u)·v(2−v)`, nought on the loop and one at
|
||||
/// the centre, the same on a spoke from either side.
|
||||
/// control points clear of the rows the conditions hold (n + 1 rows from a G^n edge), so the loop and the conditions across it do not move. One patch
|
||||
/// rises by `16·u(1−u)·v(1−v)`; patches about a centre by `u(2−u)·v(2−v)`, nought on the loop and
|
||||
/// one at the centre, the same on a spoke from either side.
|
||||
fn crowned(mut quilt: Quilt, stations: &[Station], crown: f64, size: f64) -> SplineResult<Quilt> {
|
||||
if crown == 0.0 {
|
||||
return Ok(quilt);
|
||||
}
|
||||
let held = 1 + stations.iter().map(|s| if s.across.bend.is_some() { 2 } else { 1 }).max().unwrap_or(0);
|
||||
let held = 1 + stations.iter().map(Station::order).max().unwrap_or(0);
|
||||
let about_centre = !quilt.spokes.is_empty();
|
||||
let bump = |x: f64, y: f64| if about_centre { x * (2.0 - x) * y * (2.0 - y) } else { 16.0 * x * (1.0 - x) * y * (1.0 - y) };
|
||||
let mut moved = std::collections::HashMap::new();
|
||||
@@ -983,3 +1297,7 @@ fn crowned(mut quilt: Quilt, stations: &[Station], crown: f64, size: f64) -> Spl
|
||||
}
|
||||
Ok(quilt)
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
#[path = "fill_iso_tests.rs"]
|
||||
mod tests;
|
||||
|
||||
@@ -0,0 +1,219 @@
|
||||
use super::*;
|
||||
use crate::brep::eval::SurfaceEval;
|
||||
use crate::brep::{Curve, Frame, Surface};
|
||||
use crate::linkage::surface_jet;
|
||||
|
||||
const R: f64 = 10.0;
|
||||
/// The window's half-widths in longitude and latitude.
|
||||
const ALPHA: f64 = 0.5;
|
||||
const BETA: f64 = 0.35;
|
||||
|
||||
fn spec() -> IsoSpec<'static> {
|
||||
IsoSpec { tol: 1e-6, angle_tol: 0.05f64.to_radians(), crown: 0.0, points: &[] }
|
||||
}
|
||||
|
||||
/// The arc of the circle about `frame`'s z of radius `r` from angle `a0` to `a1`, on its NURBS
|
||||
/// form (rational quadratic quarter arcs, as SMLib holds circles); run backwards when `reverse`.
|
||||
fn arc(frame: Frame, r: f64, [a0, a1]: [f64; 2], reverse: bool) -> IsoEdge {
|
||||
let (form, [t0, t1]) = crate::brep::form::curve(&Curve::Circle { frame, radius: r }, [a0, a1]).unwrap();
|
||||
let Curve::RationalSpline { knots, homogeneous } = form else { panic!("a circle's form is rational") };
|
||||
let h: Vec<HPoint> = homogeneous.iter().map(|h| HPoint { x: h[0], y: h[1], z: h[2], w: h[3] }).collect();
|
||||
let p = knots.len() - h.len() - 1;
|
||||
let (mut knots, mut h) = segment(p, &knots, &h, t0, t1).unwrap();
|
||||
if reverse {
|
||||
knots = knots.iter().rev().map(|k| t0 + t1 - k).collect();
|
||||
h.reverse();
|
||||
}
|
||||
IsoEdge { degree: p, knots, points: h.iter().map(|q| q.project().unwrap()).collect(), weights: Some(h.iter().map(|q| q.w).collect()) }
|
||||
}
|
||||
|
||||
fn frame(origin: Vec3, x: Vec3, y: Vec3) -> Frame {
|
||||
Frame { origin: [origin.x, origin.y, origin.z], x_axis: [x.x, x.y, x.z], y_axis: [y.x, y.y, y.z] }
|
||||
}
|
||||
|
||||
/// A window in a sphere of radius `R` about the origin: longitudes `±ALPHA`, latitudes `±BETA`;
|
||||
/// the bottom latitude east, the east meridian up, the top latitude west, the west meridian down.
|
||||
fn window() -> Vec<IsoEdge> {
|
||||
let radial = |u: f64| Vec3::new(u.cos(), u.sin(), 0.0);
|
||||
let latitude = |v: f64| frame(Vec3::Z * (R * v.sin()), Vec3::X, Vec3::Y);
|
||||
let meridian = |u: f64| frame(Vec3::ZERO, radial(u), Vec3::Z);
|
||||
vec![
|
||||
arc(latitude(-BETA), R * BETA.cos(), [-ALPHA, ALPHA], false),
|
||||
arc(meridian(ALPHA), R, [-BETA, BETA], false),
|
||||
arc(latitude(BETA), R * BETA.cos(), [-ALPHA, ALPHA], true),
|
||||
arc(meridian(-ALPHA), R, [-BETA, BETA], true),
|
||||
]
|
||||
}
|
||||
|
||||
/// The sphere round the window, met with continuity `order` (1 … 4): its outward normal, the way
|
||||
/// it leaves into the window, and its jet.
|
||||
fn sphere_beside(edges: &[IsoEdge], order: usize) -> impl FnMut(usize, &[f64]) -> Result<Vec<Option<Across>>, String> + '_ {
|
||||
let sphere = SurfaceEval::new(&Surface::Sphere { frame: frame(Vec3::ZERO, Vec3::X, Vec3::Y), radius: R }).unwrap();
|
||||
move |k, ts| {
|
||||
let e = &edges[k];
|
||||
ts.iter()
|
||||
.map(|t| {
|
||||
let p = eval_edge_h(e, *t).unwrap().project().unwrap();
|
||||
let uv = [p.y.atan2(p.x), (p.z / R).clamp(-1.0, 1.0).asin()];
|
||||
let (_, su, sv) = sphere.derivs(uv[0], uv[1]).unwrap();
|
||||
let normal = su.cross(sv).normalized().unwrap();
|
||||
let leave = match k {
|
||||
0 => sv.normalized().unwrap(),
|
||||
1 => su.normalized().unwrap() * -1.0,
|
||||
2 => sv.normalized().unwrap() * -1.0,
|
||||
_ => su.normalized().unwrap(),
|
||||
};
|
||||
let jet = (order >= 2).then(|| surface_jet(&sphere, uv, normal, leave, order).unwrap());
|
||||
Ok(Some(Across { leave, normal, bend: None, jet }))
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
}
|
||||
|
||||
/// How far the fill strays from the sphere, over a grid of each patch.
|
||||
fn off_sphere(fill: &IsoFill) -> f64 {
|
||||
let mut worst = 0.0f64;
|
||||
for p in &fill.patches {
|
||||
for i in 0..=16 {
|
||||
for j in 0..=16 {
|
||||
let q = p.surface.eval(i as f64 / 16.0, j as f64 / 16.0).unwrap();
|
||||
worst = worst.max((q.len() - R).abs());
|
||||
}
|
||||
}
|
||||
}
|
||||
worst
|
||||
}
|
||||
|
||||
/// The fill's one patch as a surface record, to be evaluated independently of the fill.
|
||||
fn record(fill: &IsoFill) -> SurfaceEval {
|
||||
let s = &fill.patches[0].surface;
|
||||
let [nu, nv] = s.count();
|
||||
let homogeneous = (0..nu).map(|i| (0..nv).map(|j| s.point(i, j)).map(|h| [h.x, h.y, h.z, h.w]).collect()).collect();
|
||||
let [u_knots, v_knots] = s.knots().clone();
|
||||
SurfaceEval::new(&Surface::RationalSpline { u_knots, v_knots, homogeneous }).unwrap()
|
||||
}
|
||||
|
||||
/// The largest jump of order `k` (2 … 4) across the window's edges, clear of the corners: the
|
||||
/// fill's linkage curves (its exact partials across each edge) against the sphere's jet of order
|
||||
/// `k`, per unit length of their foot, as `crate::linkage` defines it.
|
||||
fn jump(fill: &IsoFill, edges: &[IsoEdge], k: usize) -> f64 {
|
||||
let surface = record(fill);
|
||||
let mut beside = sphere_beside(edges, k);
|
||||
let mut worst = 0.0f64;
|
||||
for (side, edge) in edges.iter().enumerate() {
|
||||
let [a, b] = edge.domain();
|
||||
for i in 1..40 {
|
||||
let f = i as f64 / 40.0;
|
||||
if !(CORNER_SHARE..=1.0 - CORNER_SHARE).contains(&f) {
|
||||
continue;
|
||||
}
|
||||
let across = beside(side, &[a + (b - a) * f]).unwrap()[0].unwrap();
|
||||
let (uv, (du, dv, sign)) = station_on(side, f);
|
||||
let d = surface.partials(uv[0], uv[1]).unwrap();
|
||||
let c: Vec<Vec3> = (1..=k).map(|m| d[m * du][m * dv] * sign.powi(m as i32)).collect();
|
||||
let (n, t, w) = (across.normal, across.normal.cross(across.leave), across.leave);
|
||||
let x: Vec<[f64; 2]> = c.iter().map(|d| [d.dot(t), d.dot(w)]).collect();
|
||||
let want = across.jet.unwrap().targets(&x).unwrap()[k - 2];
|
||||
let foot = (c[0] - n * c[0].dot(n)).len();
|
||||
worst = worst.max(((n.dot(c[k - 1]) - want) / foot.powi(k as i32)).abs());
|
||||
}
|
||||
}
|
||||
worst
|
||||
}
|
||||
|
||||
/// A window in a sphere, bounded by four circular arcs (rational edges), filled G1 … G4. Each
|
||||
/// fill of order n holds its jumps of orders 2 … n, by its own report and by an independent
|
||||
/// reading of its record; the jump of order n + 1 it leaves to its energy is many times larger.
|
||||
/// (The sphere itself is not in the fill's space, whose weights are a Coons blend and the
|
||||
/// sphere's a product, so the orders a fill is not asked for do leave jumps.)
|
||||
#[test]
|
||||
fn a_sphere_window_on_rational_arcs_holds_each_order_asked_and_not_the_next() {
|
||||
let edges = window();
|
||||
for order in 1..=4 {
|
||||
let fill = iso_fill(&edges, &mut sphere_beside(&edges, order), &[false; 4], &spec()).unwrap();
|
||||
let a = fill.achieved;
|
||||
let measured: Vec<f64> = (2..=4).map(|k| jump(&fill, &edges, k)).collect();
|
||||
eprintln!(
|
||||
"G{order}: degree {}, {} control points a side, normals within {:.2e}°, reported jumps {:?}, read off the record (G2, G3, G4) {measured:?}, off the sphere {:.2e}",
|
||||
fill.patches[0].surface.degree()[0],
|
||||
fill.patches[0].surface.count()[0],
|
||||
a.normal_deviation.to_degrees(),
|
||||
a.jumps,
|
||||
off_sphere(&fill)
|
||||
);
|
||||
assert!(a.normal_deviation <= spec().angle_tol);
|
||||
// Each order held within the angle tolerance over the window's size, by the fill's own
|
||||
// report and by the independent reading of its record (at other stations, so within
|
||||
// twice that).
|
||||
let size = {
|
||||
let mut all = edges.iter().flat_map(|e| e.points.iter().copied());
|
||||
let first = all.next().unwrap();
|
||||
let (lo, hi) = all.fold((first, first), |(lo, hi), p| (lo.min(p), hi.max(p)));
|
||||
lo.dist(hi)
|
||||
};
|
||||
for k in 2..=order {
|
||||
let (reported, read) = (turn(a.jumps[k], k, size), turn(measured[k - 2], k, size));
|
||||
assert!(reported <= spec().angle_tol, "G{order}: order {k} reported {reported:e}");
|
||||
assert!(read <= 2.0 * spec().angle_tol, "G{order}: order {k} read {read:e}");
|
||||
}
|
||||
if (2..4).contains(&order) {
|
||||
let (held, next) = (measured[order - 2], measured[order - 1]);
|
||||
assert!(next > 100.0 * held.max(1e-12), "G{order} leaves order {} to the energy: {next:e} vs {held:e}", order + 1);
|
||||
}
|
||||
// Not the sphere (nothing asks for it), but a cap over the window, not a dent.
|
||||
assert!(off_sphere(&fill) < 0.05 * R, "G{order} near the sphere: {}", off_sphere(&fill));
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn hostile_edges_are_refused() {
|
||||
let edges = window();
|
||||
let ok = |e: &[IsoEdge]| iso_fill(e, &mut sphere_beside(e, 2), &[false; 4], &spec());
|
||||
let mut bad = edges.clone();
|
||||
bad[1].weights = Some(vec![-1.0; bad[1].points.len()]);
|
||||
assert!(ok(&bad).is_err(), "a negative weight");
|
||||
let mut bad = edges.clone();
|
||||
bad[1].weights = Some(vec![1.0]);
|
||||
assert!(ok(&bad).is_err(), "a weight count that does not match");
|
||||
let mut bad = edges.clone();
|
||||
if let Some(w) = bad[0].weights.as_mut() {
|
||||
w[0] = f64::NAN;
|
||||
}
|
||||
assert!(ok(&bad).is_err(), "a weight that is not a number");
|
||||
// One edge's start weight doubled alone (its curve changes, and its corner disagrees).
|
||||
let mut bad = edges.clone();
|
||||
if let Some(w) = bad[0].weights.as_mut() {
|
||||
w[0] *= 2.0;
|
||||
}
|
||||
assert!(matches!(ok(&bad), Err(SplineError::Unsupported(_))), "corner weights that cannot agree");
|
||||
// A jet above G4.
|
||||
let mut beside = sphere_beside(&edges, 4);
|
||||
let mut too_high = |k: usize, ts: &[f64]| -> Result<Vec<Option<Across>>, String> {
|
||||
let mut got = beside(k, ts)?;
|
||||
for a in got.iter_mut().flatten() {
|
||||
if let Some(j) = a.jet.as_mut() {
|
||||
j.order = 5;
|
||||
}
|
||||
}
|
||||
Ok(got)
|
||||
};
|
||||
assert!(iso_fill(&edges, &mut too_high, &[false; 4], &spec()).is_err());
|
||||
}
|
||||
|
||||
/// The rational basis's derivatives (the Leibniz rule on `w N N = R W`) against central
|
||||
/// differences of the derivative one order lower, on the sphere window's patch.
|
||||
#[test]
|
||||
fn a_rational_patch_s_derivatives_are_its_differences() {
|
||||
let edges = window();
|
||||
let q = Quilt::new(&edges, &[], 3).unwrap();
|
||||
let p = &q.patches[0];
|
||||
assert!(p.rational);
|
||||
let (u, v, h) = (0.3, 0.6, 1e-5);
|
||||
for (du, dv) in [(1, 0), (0, 1), (2, 0), (0, 2), (1, 1), (3, 0), (0, 3)] {
|
||||
let exact = p.eval(u, v, du, dv).unwrap();
|
||||
let (lu, lv) = if du > 0 { (du - 1, dv) } else { (du, dv - 1) };
|
||||
let (a, b) = if du > 0 { ((u + h, v), (u - h, v)) } else { ((u, v + h), (u, v - h)) };
|
||||
let differences = (p.eval(a.0, a.1, lu, lv).unwrap() - p.eval(b.0, b.1, lu, lv).unwrap()) * (0.5 / h);
|
||||
assert!(exact.dist(differences) < 1e-5 * exact.len().max(1.0), "∂({du}, {dv}): {exact:?} vs {differences:?}");
|
||||
}
|
||||
}
|
||||
@@ -36,7 +36,13 @@ fn setback_hole(r: f64, s: f64, g2: bool) -> impl FnMut(usize, &[f64]) -> Result
|
||||
_ => return Err(format!("no side {k}")),
|
||||
};
|
||||
let up = Vec3::new(1.0, 1.0, 1.0).normalized().unwrap_or(Vec3::Z);
|
||||
Ok(Station { p, up, across: Some(Across { leave, normal, bend: g2.then_some(Vec3::ZERO) }), corner: *t == 0.0, clash: false })
|
||||
Ok(Station {
|
||||
p,
|
||||
up,
|
||||
across: Some(Across { leave, normal, bend: g2.then_some(Vec3::ZERO), jet: None }),
|
||||
corner: *t == 0.0,
|
||||
clash: false,
|
||||
})
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
@@ -110,7 +116,7 @@ fn a_hole_in_a_cylinder_is_filled_with_the_cylinder_to_the_second_order() {
|
||||
// over the radius, bending away from the normal.
|
||||
let across = leave - Vec3::X * leave.x;
|
||||
let bend = normal * (-across.len2() / r);
|
||||
Ok(Station { p, up: normal, across: Some(Across { leave, normal, bend: Some(bend) }), corner: false, clash: false })
|
||||
Ok(Station { p, up: normal, across: Some(Across { leave, normal, bend: Some(bend), jet: None }), corner: false, clash: false })
|
||||
})
|
||||
.collect()
|
||||
};
|
||||
@@ -131,7 +137,7 @@ fn hostile_matched_fills_are_refused() {
|
||||
.map(|t| {
|
||||
let a = t * std::f64::consts::TAU / 3.0;
|
||||
let p = Vec3::new(a.cos(), a.sin(), 0.0);
|
||||
Station { p, up: Vec3::Z, across: Some(Across { leave: Vec3::Z, normal: p, bend: None }), corner: false, clash: false }
|
||||
Station { p, up: Vec3::Z, across: Some(Across { leave: Vec3::Z, normal: p, bend: None, jet: None }), corner: false, clash: false }
|
||||
})
|
||||
.collect())
|
||||
};
|
||||
|
||||
@@ -24,6 +24,7 @@ pub mod fill_general;
|
||||
pub mod fill_iso;
|
||||
pub mod fill_matched;
|
||||
mod grid;
|
||||
pub mod linkage;
|
||||
mod mesh;
|
||||
pub mod mesh_distance;
|
||||
pub mod network;
|
||||
|
||||
@@ -0,0 +1,334 @@
|
||||
//! Geometric continuity of order `n` (G0 … G4) between two surfaces that share a curve, in the
|
||||
//! linkage-curve formulation, as the fill imposes it and as a seam is measured.
|
||||
//!
|
||||
//! **The criterion.** Two surfaces meeting along a curve are G^n there exactly when the
|
||||
//! *linkage curves* crossing it — a curve on one surface continued by a curve on the other — are
|
||||
//! G^n curves (Pegna & Wolter, "Geometrical criteria to guarantee curvature continuity of blend
|
||||
//! surfaces", J. Mech. Des. 114, 1992, for n = 2; Hermann, Lukács & Wolter, "Geometrical
|
||||
//! criteria on the higher order smoothness of composite surfaces", CAGD 16, 1999, for any n).
|
||||
//! In a form that needs no reparameterisation: write the neighbour, near a point `p` of the
|
||||
//! curve, as a height `H` over its tangent plane along its unit normal `n`,
|
||||
//! `q = p + x + H(x)·n` for `x` in the tangent plane, so that `H(0) = 0` and `∇H(0) = 0`. A
|
||||
//! curve `c(σ)` from `p` lies on a surface meeting the neighbour with G^n exactly when its height
|
||||
//! above the neighbour, `n·(c(σ) − p) − H(x(σ))` with `x(σ)` the part of `c(σ) − p` in the
|
||||
//! tangent plane, vanishes to order `n` in σ. With `H(x) = Σ_{m ≥ 2} H_m(x) / m!` (`H_m` the
|
||||
//! m-th derivative of `H` at `p`, a form of degree m) and `x(σ) = Σ_j x_j σʲ / j!`
|
||||
//! (`x_j` the tangent-plane part of `c⁽ʲ⁾`), the condition of order `k` is
|
||||
//!
|
||||
//! `n·c⁽ᵏ⁾ = k! · [σᵏ] Σ_{m=2}^{k} H_m(x(σ)) / m!`,
|
||||
//!
|
||||
//! whose right-hand side reads `c′ … c⁽ᵏ⁻¹⁾` only: order 1 is `n·c′ = 0` (the tangent plane),
|
||||
//! order 2 `n·c″ = H₂(x₁)` (the normal curvature in the curve's direction), order 3
|
||||
//! `n·c‴ = 3·H₂(x₁, x₂) + H₃(x₁)` (the rate of change of that curvature), and order 4 its rate of
|
||||
//! change in turn. For a unit `a`, `H₂(a)` is the curvature of the normal section along `a` and
|
||||
//! `H₃(a)`, `H₄(a)` are its first and second derivatives by arc length.
|
||||
//!
|
||||
//! **The jet.** [`surface_jet`] reads `H₂ … H_n` off a surface record's exact partial
|
||||
//! derivatives: along a tangent vector `a`, the curve on the surface whose foot in the tangent
|
||||
//! plane runs straight along `a` at unit rate has `n·γ⁽ᵐ⁾(0) = H_m(a)`; its parameter series
|
||||
//! `(u(τ), v(τ))` is found order by order (the part of each coefficient of `S(u(τ), v(τ))` in the
|
||||
//! tangent plane is `a`, then zero), and `H_m` as a binary form from its values along
|
||||
//! `t + λ·w` at `m + 1` values of λ (a polynomial of degree m in λ, so exactly).
|
||||
//!
|
||||
//! Input is hostile: orders are capped at [`MAX_CONTINUITY`], and every evaluation is checked.
|
||||
|
||||
use crate::Vec3;
|
||||
use crate::brep::eval::SurfaceEval;
|
||||
use crate::tess::partials::{MAX_PARTIAL_ORDER, Partials};
|
||||
|
||||
/// The highest continuity order handled: G4, the highest order of the surfaces' exact partials.
|
||||
pub const MAX_CONTINUITY: usize = MAX_PARTIAL_ORDER;
|
||||
|
||||
/// A surface near a point as a height over its tangent plane, to order `order` (2 … 4): the
|
||||
/// forms `H_m` for `m = 2 ..= order` on the frame `t = n × w`, `w` (both unit, in the tangent
|
||||
/// plane; `w` the direction the surface leaves the shared curve in). `forms[m − 2][j]` is the
|
||||
/// coefficient of `x_t^(m−j)·x_w^j` in `H_m(x)`.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub struct Jet {
|
||||
pub order: usize,
|
||||
pub forms: [[f64; MAX_CONTINUITY + 1]; MAX_CONTINUITY - 1],
|
||||
}
|
||||
|
||||
impl Jet {
|
||||
/// `H_m(a)` for a tangent vector given by its components `(a_t, a_w)`.
|
||||
pub fn form(&self, m: usize, a: [f64; 2]) -> Option<f64> {
|
||||
let coefficients = self.forms.get(m.checked_sub(2)?).filter(|_| (2..=self.order).contains(&m))?;
|
||||
Some((0..=m).map(|j| coefficients.get(j).copied().unwrap_or(0.0) * a[0].powi((m - j) as i32) * a[1].powi(j as i32)).sum())
|
||||
}
|
||||
|
||||
/// Whether every coefficient is finite.
|
||||
pub fn is_finite(&self) -> bool {
|
||||
self.forms.iter().flatten().all(|x| x.is_finite())
|
||||
}
|
||||
|
||||
/// What orders `2 ..= self.order` ask of `n·c⁽ᵏ⁾` (module docs), for a curve whose first
|
||||
/// derivatives have the tangent-plane components `x[j − 1] = (x_j·t, x_j·w)`: element
|
||||
/// `k − 2` is order k's (it reads `x[0 .. k − 1]`). `None` when too few derivatives are given.
|
||||
pub fn targets(&self, x: &[[f64; 2]]) -> Option<Vec<f64>> {
|
||||
Some(self.linearised(x)?.into_iter().map(|t| t.value).collect())
|
||||
}
|
||||
|
||||
/// [`Jet::targets`] with their exact gradients: element `k − 2` is order k's target and its
|
||||
/// derivative with respect to each component of each `x_j` (`gradient[j − 1]`), so that the
|
||||
/// condition `n·c⁽ᵏ⁾ = target(x)` can be linearised about a curve (Newton's method).
|
||||
pub fn linearised(&self, x: &[[f64; 2]]) -> Option<Vec<Target>> {
|
||||
let top = self.order;
|
||||
if x.len() + 1 < top || top > MAX_CONTINUITY {
|
||||
return None;
|
||||
}
|
||||
// x(σ) as truncated Taylor series in each component (no constant term), each coefficient
|
||||
// a dual number seeded with its own component.
|
||||
let mut series = [vec![Dual::ZERO; top + 1], vec![Dual::ZERO; top + 1]];
|
||||
for (j, xj) in x.iter().enumerate().take(top.saturating_sub(1)) {
|
||||
for (c, (s, value)) in series.iter_mut().zip(xj).enumerate() {
|
||||
if let Some(slot) = s.get_mut(j + 1) {
|
||||
*slot = Dual::seed(*value, 2 * j + c).scale(1.0 / factorial(j + 1));
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut height = vec![Dual::ZERO; top + 1];
|
||||
for m in 2..=top {
|
||||
let coefficients = self.forms.get(m - 2)?;
|
||||
for (j, c) in coefficients.iter().enumerate().take(m + 1) {
|
||||
let monomial = mul(&power(&series[0], m - j, top), &power(&series[1], j, top), top);
|
||||
for (h, x) in height.iter_mut().zip(&monomial) {
|
||||
*h = *h + x.scale(c / factorial(m));
|
||||
}
|
||||
}
|
||||
}
|
||||
Some(
|
||||
(2..=top)
|
||||
.map(|k| {
|
||||
let d = height.get(k).copied().unwrap_or(Dual::ZERO).scale(factorial(k));
|
||||
let mut gradient = [[0.0; 2]; MAX_CONTINUITY - 1];
|
||||
for (j, g) in gradient.iter_mut().enumerate() {
|
||||
*g = [d.grad[2 * j], d.grad[2 * j + 1]];
|
||||
}
|
||||
Target { value: d.value, gradient }
|
||||
})
|
||||
.collect(),
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
/// One order's target and its gradient (see [`Jet::linearised`]).
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub struct Target {
|
||||
pub value: f64,
|
||||
/// `gradient[j − 1]`: the derivative with respect to `(x_j·t, x_j·w)`.
|
||||
pub gradient: [[f64; 2]; MAX_CONTINUITY - 1],
|
||||
}
|
||||
|
||||
/// Gradients carried: two components of each of `x₁ … x₃`.
|
||||
const SEEDS: usize = 2 * (MAX_CONTINUITY - 1);
|
||||
|
||||
/// A number with its gradient with respect to the [`SEEDS`] components (forward differentiation:
|
||||
/// the targets are polynomials in them, so the gradient is exact).
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
struct Dual {
|
||||
value: f64,
|
||||
grad: [f64; SEEDS],
|
||||
}
|
||||
|
||||
impl Dual {
|
||||
const ZERO: Dual = Dual { value: 0.0, grad: [0.0; SEEDS] };
|
||||
const ONE: Dual = Dual { value: 1.0, grad: [0.0; SEEDS] };
|
||||
|
||||
fn seed(value: f64, at: usize) -> Dual {
|
||||
let mut grad = [0.0; SEEDS];
|
||||
if let Some(g) = grad.get_mut(at) {
|
||||
*g = 1.0;
|
||||
}
|
||||
Dual { value, grad }
|
||||
}
|
||||
|
||||
fn scale(self, k: f64) -> Dual {
|
||||
Dual { value: self.value * k, grad: self.grad.map(|g| g * k) }
|
||||
}
|
||||
}
|
||||
|
||||
impl std::ops::Add for Dual {
|
||||
type Output = Dual;
|
||||
fn add(self, o: Dual) -> Dual {
|
||||
let mut grad = self.grad;
|
||||
for (g, h) in grad.iter_mut().zip(o.grad) {
|
||||
*g += h;
|
||||
}
|
||||
Dual { value: self.value + o.value, grad }
|
||||
}
|
||||
}
|
||||
|
||||
impl std::ops::Mul for Dual {
|
||||
type Output = Dual;
|
||||
fn mul(self, o: Dual) -> Dual {
|
||||
let mut grad = [0.0; SEEDS];
|
||||
for ((g, a), b) in grad.iter_mut().zip(self.grad).zip(o.grad) {
|
||||
*g = a * o.value + self.value * b;
|
||||
}
|
||||
Dual { value: self.value * o.value, grad }
|
||||
}
|
||||
}
|
||||
|
||||
/// `k!` as a float (small k).
|
||||
fn factorial(k: usize) -> f64 {
|
||||
(1..=k).map(|i| i as f64).product()
|
||||
}
|
||||
|
||||
/// The numbers the series below are over: plain floats and [`Dual`]s.
|
||||
trait Coefficient: Copy + std::ops::Add<Output = Self> + std::ops::Mul<Output = Self> {
|
||||
const ZERO: Self;
|
||||
const ONE: Self;
|
||||
}
|
||||
|
||||
impl Coefficient for f64 {
|
||||
const ZERO: f64 = 0.0;
|
||||
const ONE: f64 = 1.0;
|
||||
}
|
||||
|
||||
impl Coefficient for Dual {
|
||||
const ZERO: Dual = Dual::ZERO;
|
||||
const ONE: Dual = Dual::ONE;
|
||||
}
|
||||
|
||||
/// The product of two truncated series (coefficients of `σ⁰ … σ^top`).
|
||||
fn mul<T: Coefficient>(a: &[T], b: &[T], top: usize) -> Vec<T> {
|
||||
let mut out = vec![T::ZERO; top + 1];
|
||||
for (i, x) in a.iter().enumerate().take(top + 1) {
|
||||
for (j, y) in b.iter().enumerate().take(top + 1 - i) {
|
||||
if let Some(slot) = out.get_mut(i + j) {
|
||||
*slot = *slot + *x * *y;
|
||||
}
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// A truncated series to the power `e`.
|
||||
fn power<T: Coefficient>(a: &[T], e: usize, top: usize) -> Vec<T> {
|
||||
let mut one = vec![T::ZERO; top + 1];
|
||||
if let Some(c) = one.first_mut() {
|
||||
*c = T::ONE;
|
||||
}
|
||||
(0..e).fold(one, |acc, _| mul(&acc, a, top))
|
||||
}
|
||||
|
||||
/// The surface `s` near `uv` as a height along its unit normal `n` (either sign) over its tangent
|
||||
/// plane, to order `order` (2 ..= [`MAX_CONTINUITY`]), on the frame `t = n × w`, `w` (`w` a unit
|
||||
/// tangent vector): see [`Jet`] and the module docs. `None` where the surface has no tangent
|
||||
/// plane (a pole) or cannot be evaluated.
|
||||
pub fn surface_jet(s: &SurfaceEval, uv: [f64; 2], n: Vec3, w: Vec3, order: usize) -> Option<Jet> {
|
||||
if !(2..=MAX_CONTINUITY).contains(&order) {
|
||||
return None;
|
||||
}
|
||||
let d = s.partials(uv[0], uv[1])?;
|
||||
let t = n.cross(w);
|
||||
let mut forms = [[0.0; MAX_CONTINUITY + 1]; MAX_CONTINUITY - 1];
|
||||
// H_m(t + λ w) is a polynomial of degree m in λ: its coefficients are the form's.
|
||||
const NODES: [f64; MAX_CONTINUITY + 1] = [0.0, 1.0, -1.0, 2.0, -2.0];
|
||||
let along: Vec<Vec<f64>> = NODES.iter().map(|l| heights(&d, n, t + w * *l, order)).collect::<Option<_>>()?;
|
||||
for m in 2..=order {
|
||||
let values: Vec<f64> = along.iter().take(m + 1).map(|h| h.get(m).copied().unwrap_or(f64::NAN)).collect();
|
||||
let coefficients = polynomial_through(NODES.get(..=m)?, &values)?;
|
||||
let slot = forms.get_mut(m - 2)?;
|
||||
for (x, c) in slot.iter_mut().zip(coefficients) {
|
||||
*x = c;
|
||||
}
|
||||
}
|
||||
let jet = Jet { order, forms };
|
||||
jet.is_finite().then_some(jet)
|
||||
}
|
||||
|
||||
/// `H_m(a)` for `m = 0 ..= order` from the partials `d` (`d[i][j] = S^(i,j)`): the curve on
|
||||
/// the surface whose foot runs straight along `a` (module docs), its parameter series solved
|
||||
/// order by order. `None` where the tangent plane degenerates.
|
||||
fn heights(d: &Partials, n: Vec3, a: Vec3, order: usize) -> Option<Vec<f64>> {
|
||||
let at = |i: usize, j: usize| d.get(i).and_then(|r| r.get(j)).copied();
|
||||
let (su, sv) = (at(1, 0)?, at(0, 1)?);
|
||||
let (e, f, g) = (su.dot(su), su.dot(sv), sv.dot(sv));
|
||||
let det = e * g - f * f;
|
||||
// A pole: one parameter direction collapses against the other (or the two run together).
|
||||
if !(e.min(g) > 1e-20 * e.max(g) && det > 1e-12 * e * g) {
|
||||
return None;
|
||||
}
|
||||
// The rates (du, dv) whose image S_u du + S_v dv has tangent-plane part `b`.
|
||||
let rates = |b: Vec3| {
|
||||
let (bu, bv) = (b.dot(su), b.dot(sv));
|
||||
((g * bu - f * bv) / det, (e * bv - f * bu) / det)
|
||||
};
|
||||
let tangential = |x: Vec3| x - n * x.dot(n);
|
||||
let mut du = vec![0.0; order + 1];
|
||||
let mut dv = vec![0.0; order + 1];
|
||||
for k in 1..=order {
|
||||
// The coefficient of τᵏ with this order's rates still nought, then the rates that make
|
||||
// its tangent-plane part `a` (k = 1) or zero.
|
||||
let have = composed(d, &du, &dv, order)?.get(k).copied()?;
|
||||
let want = if k == 1 { tangential(a) } else { Vec3::ZERO };
|
||||
let (x, y) = rates(want - tangential(have));
|
||||
*du.get_mut(k)? = x;
|
||||
*dv.get_mut(k)? = y;
|
||||
}
|
||||
let c = composed(d, &du, &dv, order)?;
|
||||
Some(c.iter().enumerate().map(|(k, x)| factorial(k) * n.dot(*x)).collect())
|
||||
}
|
||||
|
||||
/// The Taylor coefficients of `S(u₀ + δu(τ), v₀ + δv(τ)) − S(u₀, v₀)` to order `order`, for the
|
||||
/// series `δu`, `δv` (no constant term) and the partials `d`.
|
||||
fn composed(d: &Partials, du: &[f64], dv: &[f64], order: usize) -> Option<Vec<Vec3>> {
|
||||
let mut out = vec![Vec3::ZERO; order + 1];
|
||||
for (i, row) in d.iter().enumerate() {
|
||||
for (j, partial) in row.iter().enumerate() {
|
||||
if i + j == 0 || i + j > order {
|
||||
continue;
|
||||
}
|
||||
let term = mul(&power(du, i, order), &power(dv, j, order), order);
|
||||
let scale = 1.0 / (factorial(i) * factorial(j));
|
||||
for (o, x) in out.iter_mut().zip(&term) {
|
||||
*o += *partial * (x * scale);
|
||||
}
|
||||
}
|
||||
}
|
||||
out.iter().all(|p| p.is_finite()).then_some(out)
|
||||
}
|
||||
|
||||
/// The coefficients `c₀ … c_m` of the polynomial of degree `m = nodes.len() − 1` through
|
||||
/// `(nodes[i], values[i])` (Newton's divided differences, expanded).
|
||||
fn polynomial_through(nodes: &[f64], values: &[f64]) -> Option<Vec<f64>> {
|
||||
let m = nodes.len();
|
||||
if values.len() != m || m == 0 {
|
||||
return None;
|
||||
}
|
||||
let mut dd = values.to_vec();
|
||||
for level in 1..m {
|
||||
for i in (level..m).rev() {
|
||||
let span = nodes.get(i)? - nodes.get(i - level)?;
|
||||
if span == 0.0 {
|
||||
return None;
|
||||
}
|
||||
*dd.get_mut(i)? = (dd.get(i)? - dd.get(i - 1)?) / span;
|
||||
}
|
||||
}
|
||||
// Horner on the Newton form: p(λ) = dd₀ + (λ − x₀)(dd₁ + (λ − x₁)(…)).
|
||||
let mut c = vec![0.0; m];
|
||||
for i in (0..m).rev() {
|
||||
// c ← c·(λ − x_i) + dd_i.
|
||||
let x = *nodes.get(i)?;
|
||||
let mut next = vec![0.0; m];
|
||||
for (k, ck) in c.iter().enumerate() {
|
||||
if let Some(slot) = next.get_mut(k + 1) {
|
||||
*slot += ck;
|
||||
}
|
||||
if let Some(slot) = next.get_mut(k) {
|
||||
*slot -= ck * x;
|
||||
}
|
||||
}
|
||||
if let Some(slot) = next.get_mut(0) {
|
||||
*slot += dd.get(i)?;
|
||||
}
|
||||
c = next;
|
||||
}
|
||||
c.iter().all(|x| x.is_finite()).then_some(c)
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
#[path = "linkage_tests.rs"]
|
||||
mod tests;
|
||||
@@ -0,0 +1,95 @@
|
||||
use super::*;
|
||||
use crate::brep::{Frame, Surface};
|
||||
|
||||
fn frame() -> Frame {
|
||||
Frame { origin: [1.0, -2.0, 0.5], x_axis: [1.0, 0.0, 0.0], y_axis: [0.0, 1.0, 0.0] }
|
||||
}
|
||||
|
||||
/// A sphere of radius `r` as a height over its tangent plane along the outward normal:
|
||||
/// `H(x) = √(r² − |x|²) − r`, so `H₂(a) = −|a|²/r`, `H₃ = 0`, `H₄(a) = −3|a|⁴/r³`.
|
||||
#[test]
|
||||
fn a_sphere_is_a_height_whose_forms_are_its_radius_s() {
|
||||
let r = 7.0;
|
||||
let s = SurfaceEval::new(&Surface::Sphere { frame: frame(), radius: r }).unwrap();
|
||||
for uv in [[0.3, 0.2], [2.0, -0.9], [4.5, 1.1]] {
|
||||
let (_, su, sv) = s.derivs(uv[0], uv[1]).unwrap();
|
||||
let n = su.cross(sv).normalized().unwrap();
|
||||
let w = (su + sv * 0.5).normalized().unwrap();
|
||||
let jet = surface_jet(&s, uv, n, w, 4).unwrap();
|
||||
for (a, len2) in [([1.0, 0.0], 1.0), ([0.0, 1.0], 1.0), ([1.0, 1.0], 2.0), ([0.6, -0.8], 1.0)] {
|
||||
assert!((jet.form(2, a).unwrap() + len2 / r).abs() < 1e-12, "H2 at {uv:?} along {a:?}");
|
||||
assert!(jet.form(3, a).unwrap().abs() < 1e-11, "H3 at {uv:?} along {a:?}");
|
||||
assert!((jet.form(4, a).unwrap() + 3.0 * len2 * len2 / r.powi(3)).abs() < 1e-11, "H4 at {uv:?} along {a:?}");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// A cylinder bends across its axis only: the height depends on the component round it, as the
|
||||
/// circle's `√(r² − x²) − r` does; along its axis it is straight.
|
||||
#[test]
|
||||
fn a_cylinder_bends_round_its_axis_and_not_along_it() {
|
||||
let r = 4.0;
|
||||
let s = SurfaceEval::new(&Surface::Cylinder { frame: frame(), radius: r }).unwrap();
|
||||
let uv = [0.7, 3.0];
|
||||
let (_, su, sv) = s.derivs(uv[0], uv[1]).unwrap();
|
||||
let n = su.cross(sv).normalized().unwrap();
|
||||
// w along the axis: t = n × w runs round the cylinder.
|
||||
let jet = surface_jet(&s, uv, n, sv.normalized().unwrap(), 4).unwrap();
|
||||
let around = [1.0, 0.0];
|
||||
assert!((jet.form(2, around).unwrap().abs() - 1.0 / r).abs() < 1e-12);
|
||||
assert!((jet.form(4, around).unwrap().abs() - 3.0 / r.powi(3)).abs() < 1e-11);
|
||||
for m in 2..=4 {
|
||||
assert!(jet.form(m, [0.0, 1.0]).unwrap().abs() < 1e-12, "straight along the axis: H{m}");
|
||||
}
|
||||
}
|
||||
|
||||
/// The same cylinder on its NURBS form (rational quadratic arcs, as SMLib holds it) has the same
|
||||
/// jet: the height is the surface's, whatever its parameterisation.
|
||||
#[test]
|
||||
fn a_cylinder_and_its_nurbs_form_have_one_jet() {
|
||||
let label = Surface::Cylinder { frame: frame(), radius: 4.0 };
|
||||
let (form, _) = crate::brep::form::surface(&label, [[0.0, 2.0], [-1.0, 5.0]]).unwrap();
|
||||
let (a, b) = (SurfaceEval::new(&label).unwrap(), SurfaceEval::new(&form).unwrap());
|
||||
let p = a.point([0.9, 2.0]).unwrap();
|
||||
let uv = b.project(p, None, 1e-9).unwrap();
|
||||
let (_, su, sv) = a.derivs(0.9, 2.0).unwrap();
|
||||
let n = su.cross(sv).normalized().unwrap();
|
||||
let w = (su.normalized().unwrap() * 0.6 + sv * 0.8).normalized().unwrap();
|
||||
let (ja, jb) = (surface_jet(&a, [0.9, 2.0], n, w, 4).unwrap(), surface_jet(&b, uv, n, w, 4).unwrap());
|
||||
for m in 2..=4 {
|
||||
for x in [[1.0, 0.0], [0.0, 1.0], [0.3, 0.7]] {
|
||||
let (fa, fb) = (ja.form(m, x).unwrap(), jb.form(m, x).unwrap());
|
||||
assert!((fa - fb).abs() < 1e-9 * fa.abs().max(1.0), "H{m}({x:?}): {fa} vs {fb}");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// The targets are the linkage conditions: order 2 the form in the curve's direction, order 3
|
||||
/// `3·H₂(x₁, x₂) + H₃(x₁)`.
|
||||
#[test]
|
||||
fn the_targets_are_the_height_s_taylor_terms_along_the_curve() {
|
||||
let mut forms = [[0.0; MAX_CONTINUITY + 1]; MAX_CONTINUITY - 1];
|
||||
forms[0] = [2.0, 0.5, -1.0, 0.0, 0.0];
|
||||
forms[1] = [0.3, -0.2, 0.1, 0.7, 0.0];
|
||||
let jet = Jet { order: 3, forms };
|
||||
let (x1, x2) = ([0.8, 0.4], [0.1, -0.5]);
|
||||
let got = jet.targets(&[x1, x2]).unwrap();
|
||||
let h2 = |a: [f64; 2], b: [f64; 2]| 2.0 * a[0] * b[0] + 0.25 * (a[0] * b[1] + a[1] * b[0]) - a[1] * b[1];
|
||||
assert!((got[0] - h2(x1, x1)).abs() < 1e-14);
|
||||
assert!((got[1] - (3.0 * h2(x1, x2) + jet.form(3, x1).unwrap())).abs() < 1e-14);
|
||||
assert!(jet.targets(&[x1]).is_none(), "order 3 reads c″");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn hostile_orders_and_degenerate_points_are_refused() {
|
||||
let s = SurfaceEval::new(&Surface::Sphere { frame: frame(), radius: 2.0 }).unwrap();
|
||||
let (n, w) = (Vec3::Z, Vec3::X);
|
||||
assert!(surface_jet(&s, [0.0, 0.0], Vec3::X, Vec3::Y, 5).is_none());
|
||||
assert!(surface_jet(&s, [0.0, 0.0], Vec3::X, Vec3::Y, 1).is_none());
|
||||
// At the pole the latitude's parameterisation has no tangent plane.
|
||||
assert!(surface_jet(&s, [0.0, std::f64::consts::FRAC_PI_2], n, w, 2).is_none());
|
||||
assert!(surface_jet(&s, [f64::NAN, 0.0], n, w, 2).is_none());
|
||||
let jet = Jet { order: 2, forms: [[f64::NAN; MAX_CONTINUITY + 1]; MAX_CONTINUITY - 1] };
|
||||
assert!(!jet.is_finite());
|
||||
assert!(jet.form(3, [1.0, 0.0]).is_none(), "beyond its order");
|
||||
}
|
||||
@@ -76,6 +76,30 @@ impl Linear for Vec3 {
|
||||
}
|
||||
}
|
||||
|
||||
/// Homogeneous points as they are: a rational spline's control values `(x·w, y·w, z·w, w)`
|
||||
/// combine linearly, so carrying one into a larger space (degree elevation, knot insertion) is
|
||||
/// exact on them.
|
||||
impl Linear for HPoint {
|
||||
fn zero() -> Self {
|
||||
HPoint::ZERO
|
||||
}
|
||||
fn plus(self, o: Self) -> Self {
|
||||
self + o
|
||||
}
|
||||
fn times(self, s: f64) -> Self {
|
||||
self * s
|
||||
}
|
||||
fn to_h(self) -> HPoint {
|
||||
self
|
||||
}
|
||||
fn from_h(h: HPoint) -> Self {
|
||||
h
|
||||
}
|
||||
fn is_finite(self) -> bool {
|
||||
self.x.is_finite() && self.y.is_finite() && self.z.is_finite() && self.w.is_finite()
|
||||
}
|
||||
}
|
||||
|
||||
/// A clamped B-spline space: a degree and a knot vector whose first and last `degree + 1` knots
|
||||
/// are the domain's ends, interior knots repeated at most `degree` times (the space's splines
|
||||
/// are continuous).
|
||||
|
||||
+78
-50
@@ -5,11 +5,14 @@
|
||||
//! parameters are plane coordinates, so each edge's pcurve is its exact projection onto the
|
||||
//! plane.
|
||||
//!
|
||||
//! The general fill, which also meets the faces beside the loop with G1 or G2, passes through
|
||||
//! The general fill, which also meets the faces beside the loop with G1 … G4, passes through
|
||||
//! given points, and fills loops that bound no height field, is a patch whose boundary rows are
|
||||
//! the loop's own edges (`solidcraft_geom::fill_iso`): each edge's pcurve is a straight line in
|
||||
//! (u, v). It takes loops of three or four line or spline edges; the frames across each edge are
|
||||
//! read off the neighbouring face's record surface (exact second derivatives for G2).
|
||||
//! (u, v), affine in the edge's own parameter. It takes loops of three to sixteen line, spline or
|
||||
//! rational edges (a circle's arc on the NURBS form SMLib holds it as); the frames across each
|
||||
//! edge are read off the neighbouring face's record surface, with its exact derivatives to the
|
||||
//! order asked: the neighbour as a height over its tangent plane to that order
|
||||
//! (`solidcraft_geom::linkage`).
|
||||
|
||||
use solidcraft_geom::Vec3;
|
||||
use solidcraft_geom::brep::eval::SurfaceEval;
|
||||
@@ -25,7 +28,7 @@ use crate::history::{Change, EdgeOrigin, EdgeProv, FaceProv, History};
|
||||
use crate::sew::{PartEdge, PartFace, PartUse, Parts, assemble};
|
||||
use crate::{KernelError, Result};
|
||||
|
||||
/// How far (degrees) a G1 or G2 patch's normal may turn from its neighbour's along the loop.
|
||||
/// How far (degrees) a G1 … G4 patch's normal may turn from its neighbour's along the loop.
|
||||
pub const PATCH_ANGLE_TOL_DEG: f64 = 0.05;
|
||||
|
||||
/// A spline error as the patch reports it.
|
||||
@@ -39,7 +42,7 @@ fn spline_error(e: SplineError) -> KernelError {
|
||||
/// What a patch on a loop of edges must do besides filling it.
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
pub struct PatchSpec {
|
||||
/// How every edge's neighbour is met (G0, G1 or G2).
|
||||
/// How every edge's neighbour is met (G0 … G4).
|
||||
pub continuity: Continuity,
|
||||
/// Edges met otherwise: a point on the edge and its continuity.
|
||||
pub edges: Vec<(Vec3, Continuity)>,
|
||||
@@ -77,9 +80,11 @@ pub struct PatchReport {
|
||||
/// corners where the neighbours themselves disagree, and anywhere.
|
||||
pub normal_deviation_deg: f64,
|
||||
pub normal_deviation_all_deg: f64,
|
||||
/// The largest difference (1/mm) between the patch's normal curvature across a G2 edge and
|
||||
/// the neighbour's.
|
||||
pub curvature_mismatch: f64,
|
||||
/// `jumps[k]`, k = 2 … 4: the largest jump across the edges asked for order k or more, clear
|
||||
/// of the loop's corners, of the k-th derivative of the linkage curves' height above the
|
||||
/// neighbour (1/mmᵏ⁻¹): the normal curvature (k = 2), its rate of change
|
||||
/// (3) and that rate's (4) (`solidcraft_geom::fill_iso::Achieved::jumps`).
|
||||
pub jumps: [f64; solidcraft_geom::linkage::MAX_CONTINUITY + 1],
|
||||
/// The farthest point from the surface.
|
||||
pub point_gap: f64,
|
||||
}
|
||||
@@ -100,10 +105,6 @@ pub fn patch_edges_report(body: &Body, at: &[Vec3], spec: &PatchSpec) -> Result<
|
||||
if !(solidcraft_geom::fill::MIN_CROWN..=solidcraft_geom::fill::MAX_CROWN).contains(&spec.crown) {
|
||||
return Err(KernelError::Invalid("a patch's crown must be between -0.9 and 3".into()));
|
||||
}
|
||||
let g3 = |g: Continuity| g == Continuity::G3;
|
||||
if g3(spec.continuity) || spec.edges.iter().any(|(_, g)| g3(*g)) {
|
||||
return Err(KernelError::Invalid("not supported yet: a patch meets its neighbours with g0, g1 or g2".into()));
|
||||
}
|
||||
if spec.edges.len() > 10_000 || spec.edges.iter().any(|(p, _)| !p.is_finite()) {
|
||||
return Err(KernelError::Invalid("pick at most 10000 edges, each by a finite point, for their continuity".into()));
|
||||
}
|
||||
@@ -125,7 +126,7 @@ pub fn patch_edges_report(body: &Body, at: &[Vec3], spec: &PatchSpec) -> Result<
|
||||
}
|
||||
}
|
||||
if spec.crown != 0.0 && orders.iter().all(|g| *g == Continuity::G0) {
|
||||
return Err(KernelError::Invalid("not supported yet: a crown on a patch through points with no g1 or g2 edge".into()));
|
||||
return Err(KernelError::Invalid("not supported yet: a crown on a patch through points with no g1 or higher edge".into()));
|
||||
}
|
||||
general_fill(body, &record, &chain, &orders, &spec.points, spec.crown)
|
||||
}
|
||||
@@ -151,7 +152,7 @@ fn general_fill(
|
||||
let size = body.size().max(solidcraft_geom::LINEAR_TOL);
|
||||
let mut beside: Vec<Option<Beside>> = Vec::with_capacity(chain.len());
|
||||
for (&(e, forward), order) in chain.iter().zip(orders) {
|
||||
beside.push(if *order == Continuity::G0 { None } else { Some(Beside::new(record, e, forward, *order == Continuity::G2, size)?) });
|
||||
beside.push(if *order == Continuity::G0 { None } else { Some(Beside::new(record, e, forward, *order, size)?) });
|
||||
}
|
||||
// Where a G1/G2 edge's neighbour and the one before disagree at their corner.
|
||||
let n = chain.len();
|
||||
@@ -181,7 +182,7 @@ fn general_fill(
|
||||
gap: 0.0,
|
||||
normal_deviation_deg: a.normal_deviation.to_degrees(),
|
||||
normal_deviation_all_deg: a.normal_deviation_all.to_degrees(),
|
||||
curvature_mismatch: a.curvature_mismatch,
|
||||
jumps: a.jumps,
|
||||
point_gap: a.point_deviation,
|
||||
};
|
||||
Ok((made, report))
|
||||
@@ -208,7 +209,7 @@ pub(crate) fn matched_patch(body: &Body, chain: &[(usize, bool)], curvature_tol:
|
||||
let mut beside = Vec::with_capacity(chain.len());
|
||||
for &(e, forward) in chain {
|
||||
sides.push(Side::new(&record, e, forward)?);
|
||||
beside.push(Beside::new(&record, e, forward, true, size)?);
|
||||
beside.push(Beside::new(&record, e, forward, Continuity::G2, size)?);
|
||||
}
|
||||
// A vertex of the loop is a corner where the loop turns there: the edge arriving and the one
|
||||
// leaving go different ways.
|
||||
@@ -241,44 +242,44 @@ pub(crate) fn matched_patch(body: &Body, chain: &[(usize, bool)], curvature_tol:
|
||||
Ok((made, fill))
|
||||
}
|
||||
|
||||
/// A loop edge as the patch takes it: a non-rational B-spline on the edge's interval, running
|
||||
/// the loop's way, its ends on the edge's vertices. Conics and traced curves have no such form
|
||||
/// in their own parameter (an arc's angle is no rational function of a NURBS parameter), and a
|
||||
/// rational edge none the Coons sum can blend: refused by name.
|
||||
/// A loop edge as the patch takes it: a B-spline on the edge's interval, rational when the edge
|
||||
/// is (a circle or ellipse on the NURBS form SMLib holds it as, whose parameter the edge's interval
|
||||
/// and pcurves are in), running the loop's way, its ends on the edge's vertices. A traced curve
|
||||
/// has no such form in its own parameter and is refused by name.
|
||||
fn iso_edge(record: &b::Brep, e: usize, forward: bool) -> Result<IsoEdge> {
|
||||
let edge = record.edges.get(e).ok_or_else(|| KernelError::Failed("patch: edge reference".into()))?;
|
||||
let v3 = |p: &[f64; 3]| Vec3::new(p[0], p[1], p[2]);
|
||||
let [a, z] = edge.interval;
|
||||
let (degree, knots, points): (usize, Vec<f64>, Vec<Vec3>) = match &edge.curve {
|
||||
let (degree, knots, h): (usize, Vec<f64>, Vec<HPoint>) = match edge.parameters() {
|
||||
b::Curve::Line { start, end } => {
|
||||
let (s, t) = (v3(start), v3(end));
|
||||
(1, vec![a, a, z, z], vec![s + (t - s) * a, s + (t - s) * z])
|
||||
(1, vec![a, a, z, z], vec![HPoint::weighted(s + (t - s) * a, 1.0), HPoint::weighted(s + (t - s) * z, 1.0)])
|
||||
}
|
||||
b::Curve::Polyline { points } => {
|
||||
let mut knots = vec![0.0];
|
||||
knots.extend((0..points.len()).map(|i| i as f64));
|
||||
knots.push(points.len().saturating_sub(1) as f64);
|
||||
(1, knots, points.iter().map(v3).collect())
|
||||
(1, knots, points.iter().map(|p| HPoint::weighted(v3(p), 1.0)).collect())
|
||||
}
|
||||
b::Curve::Spline { knots, points } => (knots.len().saturating_sub(points.len() + 1), knots.clone(), points.iter().map(v3).collect()),
|
||||
b::Curve::RationalSpline { knots, homogeneous } if homogeneous.windows(2).all(|w| w[0][3] == w[1][3]) => (
|
||||
b::Curve::Spline { knots, points } => {
|
||||
(knots.len().saturating_sub(points.len() + 1), knots.clone(), points.iter().map(|p| HPoint::weighted(v3(p), 1.0)).collect())
|
||||
}
|
||||
b::Curve::RationalSpline { knots, homogeneous } => (
|
||||
knots.len().saturating_sub(homogeneous.len() + 1),
|
||||
knots.clone(),
|
||||
homogeneous.iter().map(|h| Vec3::new(h[0] / h[3], h[1] / h[3], h[2] / h[3])).collect(),
|
||||
homogeneous.iter().map(|h| HPoint { x: h[0], y: h[1], z: h[2], w: h[3] }).collect(),
|
||||
),
|
||||
_ => {
|
||||
return Err(KernelError::Invalid(format!(
|
||||
"not supported yet: a general patch on edge {e}, a circle, ellipse, rational or traced curve (its pcurve on the patch has no exact form in its own parameter)"
|
||||
"not supported yet: a general patch on edge {e}, a circle or ellipse with no NURBS form, or a traced curve (its pcurve on the patch has no exact form in its own parameter)"
|
||||
)));
|
||||
}
|
||||
};
|
||||
// The piece on the edge's interval, then run the loop's way.
|
||||
let h: Vec<HPoint> = points.iter().map(|p| HPoint::weighted(*p, 1.0)).collect();
|
||||
let (mut knots, h) = solidcraft_geom::bspline::segment(degree, &knots, &h, a, z).map_err(spline_error)?;
|
||||
let mut points: Vec<Vec3> = h.iter().map(|p| p.xyz()).collect();
|
||||
let (mut knots, mut h) = solidcraft_geom::bspline::segment(degree, &knots, &h, a, z).map_err(spline_error)?;
|
||||
if !forward {
|
||||
knots = knots.iter().rev().map(|k| a + z - k).collect();
|
||||
points.reverse();
|
||||
h.reverse();
|
||||
// The interval's ends exactly (the reflection's rounding), so a piece's ends name them.
|
||||
let n = knots.len();
|
||||
for (i, k) in knots.iter_mut().enumerate() {
|
||||
@@ -289,6 +290,12 @@ fn iso_edge(record: &b::Brep, e: usize, forward: bool) -> Result<IsoEdge> {
|
||||
}
|
||||
}
|
||||
}
|
||||
let rational = h.iter().any(|q| q.w != 1.0);
|
||||
let mut points = h
|
||||
.iter()
|
||||
.map(|q| q.project().filter(|_| q.w > 0.0))
|
||||
.collect::<Option<Vec<Vec3>>>()
|
||||
.ok_or_else(|| KernelError::Invalid(format!("not supported yet: a general patch on edge {e}, whose weights are not all positive")))?;
|
||||
// Its ends exactly on the loop's vertices.
|
||||
let vertex = |end: usize| record.vertices.get(edge.vertices[end] as usize).map(|v| v3(&v.point));
|
||||
let (first, last) = if forward { (vertex(0), vertex(1)) } else { (vertex(1), vertex(0)) };
|
||||
@@ -298,27 +305,28 @@ fn iso_edge(record: &b::Brep, e: usize, forward: bool) -> Result<IsoEdge> {
|
||||
*l = q;
|
||||
}
|
||||
}
|
||||
Ok(IsoEdge { degree, knots, points })
|
||||
let weights = rational.then(|| h.iter().map(|q| q.w).collect());
|
||||
Ok(IsoEdge { degree, knots, points, weights })
|
||||
}
|
||||
|
||||
/// The face beside a G1/G2 edge, read at points of the edge: its record surface, which side of
|
||||
/// it is outward, whether its boundary runs the edge the loop's way, and whether its curvature
|
||||
/// is read (G2).
|
||||
/// The face beside a G1 … G4 edge, read at points of the edge: its record surface, which side of
|
||||
/// it is outward, whether its boundary runs the edge the loop's way, and how much of it is read
|
||||
/// (`order`: its tangent plane for G1, and its jet to that order above).
|
||||
struct Beside {
|
||||
surface: SurfaceEval,
|
||||
side: f64,
|
||||
along: bool,
|
||||
edge: usize,
|
||||
forward: bool,
|
||||
curvature: bool,
|
||||
order: Continuity,
|
||||
reach: f64,
|
||||
}
|
||||
|
||||
impl Beside {
|
||||
fn new(record: &b::Brep, e: usize, forward: bool, curvature: bool, size: f64) -> Result<Beside> {
|
||||
fn new(record: &b::Brep, e: usize, forward: bool, order: Continuity, size: f64) -> Result<Beside> {
|
||||
let edge = record.edges.get(e).ok_or_else(|| KernelError::Failed("patch: edge reference".into()))?;
|
||||
let [u] = edge.uses.as_slice() else {
|
||||
return Err(KernelError::Invalid("an edge of the loop lies between two faces: g1 and g2 need the one face beside it".into()));
|
||||
return Err(KernelError::Invalid("an edge of the loop lies between two faces: g1 and up need the one face beside it".into()));
|
||||
};
|
||||
let face = record.faces.get(u.face as usize).ok_or_else(|| KernelError::Failed("patch: face reference".into()))?;
|
||||
let surface = SurfaceEval::new(face.parameters())
|
||||
@@ -345,15 +353,17 @@ impl Beside {
|
||||
along: runs_forward == forward,
|
||||
edge: e,
|
||||
forward,
|
||||
curvature,
|
||||
order,
|
||||
reach: (size * 1e-4).max(1e-6),
|
||||
})
|
||||
}
|
||||
|
||||
/// The neighbour's conditions at parameter `t` of the loop's edge (its `IsoEdge`, run the
|
||||
/// loop's way, whose parameter is the edge's reversed where the loop runs against it):
|
||||
/// its normal, the direction it leaves the edge in, away from it, and for G2 its curvature
|
||||
/// vector along that direction. `None` where the edge has no tangent there.
|
||||
/// its normal, the direction it leaves the edge in, away from it, and for G2 and up its
|
||||
/// curvature vector along that direction and its jet to the order asked
|
||||
/// (`solidcraft_geom::linkage::surface_jet`, on the frame of the boundary and the leave
|
||||
/// direction). `None` where the edge has no tangent there.
|
||||
fn across(&self, record: &b::Brep, t: f64) -> Result<Option<Across>> {
|
||||
let edge = record.edges.get(self.edge).ok_or_else(|| KernelError::Failed("patch: edge reference".into()))?;
|
||||
let [a, z] = edge.interval;
|
||||
@@ -370,11 +380,23 @@ impl Beside {
|
||||
// Away from the face: its boundary's direction × its outward normal.
|
||||
let boundary = if self.along { tangent } else { tangent * -1.0 };
|
||||
let leave = boundary.cross(normal).normalized().ok_or_else(|| KernelError::Invalid("the boundary runs along the face's normal".into()))?;
|
||||
let bend = if self.curvature { Some(normal * normal_curvature(&self.surface, uv, normal, leave)?) } else { None };
|
||||
Ok(Some(Across { leave, normal, bend }))
|
||||
let order = continuity_order(self.order);
|
||||
if order < 2 {
|
||||
return Ok(Some(Across { leave, normal, bend: None, jet: None }));
|
||||
}
|
||||
let bend = Some(normal * normal_curvature(&self.surface, uv, normal, leave)?);
|
||||
let jet = solidcraft_geom::linkage::surface_jet(&self.surface, uv, normal, leave, order).ok_or_else(|| {
|
||||
KernelError::Invalid(format!("not supported yet: g{order} to a face whose surface has no derivatives of order {order} there"))
|
||||
})?;
|
||||
Ok(Some(Across { leave, normal, bend, jet: Some(jet) }))
|
||||
}
|
||||
}
|
||||
|
||||
/// The continuity order a [`Continuity`] asks: 0 (G0) … n, the one place it is read as a number.
|
||||
fn continuity_order(g: Continuity) -> usize {
|
||||
g.conditions().saturating_sub(1)
|
||||
}
|
||||
|
||||
/// The normal curvature of a surface at `uv` along the unit tangent `w`: `n·r''` of the unit-speed
|
||||
/// curve on it leaving that way (the second fundamental form, exact second derivatives).
|
||||
fn normal_curvature(s: &SurfaceEval, uv: [f64; 2], n: Vec3, w: Vec3) -> Result<f64> {
|
||||
@@ -783,7 +805,9 @@ fn quilt_body(body: &Body, record: &b::Brep, chain: &[(usize, bool)], fill: &Iso
|
||||
};
|
||||
}
|
||||
let id = index(edges.len())?;
|
||||
edges.push(PartEdge { curve: edge.curve.clone(), interval: [lo, hi], vertices: ends });
|
||||
// The curve the interval and the pcurve's parameter are on: a conic's
|
||||
// NURBS form.
|
||||
edges.push(PartEdge { curve: edge.parameters().clone(), interval: [lo, hi], vertices: ends });
|
||||
origin.push(Some(e));
|
||||
piece_of.insert(key, id);
|
||||
id
|
||||
@@ -802,12 +826,16 @@ fn quilt_body(body: &Body, record: &b::Brep, chain: &[(usize, bool)], fill: &Iso
|
||||
let start = middle(&mut vertices, &mut middle_of, spoke, edge, 0.5 * (edge.interval[0] + edge.interval[1]))?;
|
||||
let end = centre.ok_or_else(|| fail("centre"))?;
|
||||
let id = index(edges.len())?;
|
||||
let points = curve.points.iter().map(|p| [p.x, p.y, p.z]).collect();
|
||||
edges.push(PartEdge {
|
||||
curve: b::Curve::Spline { knots: curve.knots.clone(), points },
|
||||
interval: [0.0, 1.0],
|
||||
vertices: [start, end],
|
||||
});
|
||||
let spline = match &curve.weights {
|
||||
Some(w) => b::Curve::RationalSpline {
|
||||
knots: curve.knots.clone(),
|
||||
homogeneous: curve.points.iter().zip(w).map(|(p, w)| [p.x * w, p.y * w, p.z * w, *w]).collect(),
|
||||
},
|
||||
None => {
|
||||
b::Curve::Spline { knots: curve.knots.clone(), points: curve.points.iter().map(|p| [p.x, p.y, p.z]).collect() }
|
||||
}
|
||||
};
|
||||
edges.push(PartEdge { curve: spline, interval: [0.0, 1.0], vertices: [start, end] });
|
||||
origin.push(None);
|
||||
spoke_of.insert(spoke, id);
|
||||
id
|
||||
|
||||
@@ -195,7 +195,7 @@ fn spec(g: Continuity) -> PatchSpec {
|
||||
fn report(name: &str, r: &PatchReport) {
|
||||
eprintln!(
|
||||
"{name}: fill reports gap {:e}, normals within {:.4}° ({:.4}° near clashing corners), curvature mismatch {:e}, points within {:e}",
|
||||
r.gap, r.normal_deviation_deg, r.normal_deviation_all_deg, r.curvature_mismatch, r.point_gap
|
||||
r.gap, r.normal_deviation_deg, r.normal_deviation_all_deg, r.jumps[2], r.point_gap
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Reference in New Issue
Block a user