masi/solvecraft
1
0
Fork 0
mirror of https://github.com/bherbruck/solidcraft.git synced 2026-11-03 02:25:55 +00:00
Code Issues Packages Projects Releases 1 Wiki Activity

brep, step: the pcurve fit is gone

Browse Source 
`On::fit` and everything it used (samples, least squares, breaks, filled poles) are deleted, and
the STEP reader no longer retries a face at ten times its tolerance. A pcurve is exact (closed
form, an affine image, or a straight line along a form's parameters) or the face is refused:
"not supported yet: no exact pcurve of <curve> on <surface>".

The rig and ladder are as with the fit (333/372, 101/116; none lost; 3806, already failing, now
refused on import). 16 unit tests that build sphere meridians, torus seams, closed spline
bands, intersection curves and surfaces of revolution on their labels now refuse: those paths
need meridians labelled on their face's meridian plane and revolutions on the quarter grid.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01BnQawySUFx57iaJZq7z1Fm
This commit is contained in:
jyviko
2026-10-23 18:20:36 -04:00
parent 819928205a
commit 6cc0401fa7
2 changed files with 47 additions and 336 deletions
Download Patch File Download Diff File Expand all files Collapse all files
crates/geom/src/brep/pcurve.rs
+26 -292
View File
@@ -57,12 +57,6 @@ pub struct EdgeUse<'a> {
/// Knot spans a fitted pcurve may have (the schema takes 4096 control points).
const MAX_SPANS: usize = 4000;
/// Refinement rounds of one fit.
const MAX_ROUNDS: usize = 40;
/// Samples per knot span when a fit starts.
const PER_SPAN: usize = 4;
/// Spans a fit starts from (more where the curve has its own knots).
const START_SPANS: usize = 4;
/// Edge uses of one face.
const MAX_USES: usize = 100_000;
/// Lift checks of an exact pcurve: evenly spaced, plus one inside each of its spans up to this.
@@ -170,6 +164,19 @@ fn cut(c: Curve2, [t0, t1]: [f64; 2]) -> Result<Curve2> {
})
}
/// What a curve record is, for a message.
fn curve_kind(c: &Curve) -> &'static str {
match c {
Curve::Line { .. } => "a line",
Curve::Circle { .. } => "a circle",
Curve::Ellipse { .. } => "an ellipse",
Curve::Polyline { .. } => "a polyline",
Curve::Spline { .. } => "a B-spline",
Curve::RationalSpline { .. } => "a NURBS",
_ => "a derived curve",
}
}
/// A face's surface: its evaluator and parameter box, and the grid a point's parameters on it
/// are searched from, made once for all the face's points.
struct On {
@@ -180,18 +187,6 @@ struct On {
face: Option<[[f64; 2]; 2]>,
}
/// A curve sample: its parameter, point, surface parameters (periodic ones unwrapped along the
/// curve), whether the surface is degenerate along each direction there (a pole), and the
/// curve's point halfway to the next sample once measured.
#[derive(Clone, Copy)]
struct Sample {
t: f64,
p: Vec3,
uv: [f64; 2],
pole: [bool; 2],
mid: Option<Vec3>,
}
impl On {
fn new(s: &Surface) -> Result<Self> {
let surface = SurfaceEval::new(s)?;
@@ -199,6 +194,18 @@ impl On {
Ok(On { surface, starts: uv::Starts::default(), dom, face: None })
}
/// What the surface is, for a message.
fn surface_kind(&self) -> &'static str {
match &self.surface {
SurfaceEval::Eval(crate::tess::EvalSurface::Plane { .. }, _) => "a plane",
SurfaceEval::Cylinder { .. } => "a cylinder",
SurfaceEval::Cone { .. } => "a cone",
SurfaceEval::Sphere { .. } => "a sphere",
SurfaceEval::Torus { .. } => "a torus",
_ => "a spline surface",
}
}
fn within(self, face: Option<[[f64; 2]; 2]>) -> Self {
let face = face.filter(|f| f.iter().flatten().all(|x| x.is_finite()));
On { face, ..self }
@@ -210,7 +217,7 @@ impl On {
{
return Ok(pc);
}
self.fit(record, curve, range, tol)
Err(Error::Failed(format!("not supported yet: no exact pcurve of {} on {}", curve_kind(record), self.surface_kind())))
}
/// The closed-form pcurve where there is one (unchecked).
@@ -355,12 +362,6 @@ impl On {
Some(q)
}
fn sample(&self, curve: &CurveEval, t: f64, near: Option<[f64; 2]>, accept: f64, reach: f64) -> Result<Sample> {
let p = curve.point(t).ok_or_else(|| Error::Failed(format!("the curve has no point at {t}")))?;
let uv = self.project(p, near, accept).ok_or_else(|| self.off(p, t, near))?;
Ok(Sample { t, p, uv, pole: self.pole(uv, reach), mid: None })
}
/// Whether `uv` is at a pole along each direction, or within `reach` of one: moving that
/// parameter by half its period (or across its range) moves the point by no more than `reach`
/// at its speed there. There the parameter is as free as at the pole itself.
@@ -380,99 +381,6 @@ impl On {
out
}
/// The error for a curve point `p` (at `t`) farther than the fit allows from the surface.
fn off(&self, p: Vec3, t: f64, near: Option<[f64; 2]>) -> Error {
let d = self.project(p, near, f64::INFINITY).and_then(|q| self.surface.point(q)).map(|s| s.dist(p));
Error::Failed(format!("the curve at {t} lies {} mm off the surface", d.unwrap_or(f64::NAN)))
}
/// A cubic B-spline in (u, v) on the curve's own parameter, through the parameters of the
/// curve's ends, least squares through samples between; spans whose lift misses the curve
/// by more than half the tolerance (at the samples and halfway between them) are halved and
/// sampled twice as densely until none does.
fn fit(&self, record: &Curve, curve: &CurveEval, [t0, t1]: [f64; 2], tol: f64) -> Result<Curve2> {
let (target, accept) = (0.5 * tol, 0.5 * tol);
// Samples this near a pole take their free parameter from a neighbour, which moves
// their lift by at most this (well inside the target).
let reach = 0.25 * tol;
// Interior knots: where the curve's own parameterisation is less smooth (a pcurve is no
// smoother than its curve), else even.
let mut inner = breaks(record, Some(curve), t0, t1);
if inner.len() + 1 > MAX_SPANS / 4 {
return Err(Error::TooLarge(format!("a pcurve of a curve with {} knots", inner.len())));
}
if inner.is_empty() {
inner = (1..START_SPANS).map(|k| (t0 + (t1 - t0) * k as f64 / START_SPANS as f64, 1)).collect();
}
let bounds =
|inner: &[(f64, usize)]| -> Vec<f64> { std::iter::once(t0).chain(inner.iter().map(|k| k.0)).chain(std::iter::once(t1)).collect() };
let mut samples: Vec<Sample> = Vec::new();
let mut near = None;
for w in bounds(&inner).windows(2) {
let (a, b) = (w.first().copied().unwrap_or(t0), w.get(1).copied().unwrap_or(t1));
for k in usize::from(!samples.is_empty())..=PER_SPAN {
let t = if k == PER_SPAN { b } else { a + (b - a) * k as f64 / PER_SPAN as f64 };
let s = self.sample(curve, t, near, accept, reach)?;
near = Some(s.uv);
samples.push(s);
}
}
for _ in 0..MAX_ROUNDS {
let knots: Vec<f64> = [t0; 4].into_iter().chain(inner.iter().flat_map(|&(k, m)| std::iter::repeat_n(k, m))).chain([t1; 4]).collect();
let points = least_squares(&knots, &fill_poles(&samples)).ok_or_else(|| Error::Failed("the pcurve fit is singular".into()))?;
let pc = Curve2::Spline { knots, points };
let miss = |t: f64, p: Vec3| pc.eval(t).and_then(|q| self.surface.point(q)).map_or(f64::INFINITY, |s| s.dist(p));
let spans = inner.len() + 1;
let span_of = |t: f64| inner.partition_point(|k| k.0 <= t).min(spans - 1);
let mut bad = vec![false; spans];
for i in 0..samples.len() {
let Some(s) = samples.get(i).copied() else { break };
let mid_t = samples.get(i + 1).map(|n| 0.5 * (s.t + n.t));
let mid = s.mid.or_else(|| mid_t.and_then(|t| curve.point(t)));
if let Some(slot) = samples.get_mut(i) {
slot.mid = mid;
}
let between = match (mid_t, mid) {
(Some(t), Some(m)) => miss(t, m),
(Some(_), None) => f64::INFINITY,
_ => 0.0,
};
if !(miss(s.t, s.p) <= target && between <= target)
&& let Some(flag) = bad.get_mut(span_of(mid_t.unwrap_or(s.t)))
{
*flag = true;
}
}
if !bad.contains(&true) {
return Ok(pc);
}
// Halve the failing spans and sample them twice as densely.
let split: Vec<(f64, f64)> =
bounds(&inner).windows(2).zip(&bad).filter(|(_, b)| **b).filter_map(|(w, _)| Some((*w.first()?, *w.get(1)?))).collect();
if inner.iter().map(|k| k.1).sum::<usize>() + split.len() + 1 > MAX_SPANS {
break;
}
inner.extend(split.iter().map(|(a, b)| (0.5 * (a + b), 1)));
inner.sort_by(|a, b| a.0.total_cmp(&b.0));
let inside = |t: f64| split.iter().any(|&(a, b)| t >= a && t < b);
let mut dense = Vec::with_capacity(samples.len() * 2);
for i in 0..samples.len() {
let Some(s) = samples.get(i).copied() else { break };
match (inside(s.t), samples.get(i + 1), s.mid) {
(true, Some(n), Some(m)) => {
let t = 0.5 * (s.t + n.t);
let uv = self.project(m, Some(s.uv), accept).ok_or_else(|| self.off(m, t, Some(s.uv)))?;
dense.push(Sample { mid: None, ..s });
dense.push(Sample { t, p: m, uv, pole: self.pole(uv, reach), mid: None });
}
_ => dense.push(s),
}
}
samples = dense;
}
Err(Error::Failed(format!("no pcurve within {tol} mm in {MAX_SPANS} spans")))
}
/// The loop's uses moved by whole periods so each starts where the one before it ends (in
/// loop order and direction), then all together into the surface's base period where the
/// loop's lowest point allows. Where uses meet at a pole (within `tol`) the periodic
@@ -570,59 +478,6 @@ impl On {
}
}
/// The interior parameters where `curve` is less than twice continuously differentiable, with
/// the knot multiplicity a cubic needs there: a spline's knots by their multiplicity, and the
/// vertices of a polyline or an intersection curve's leader polyline (its parameter turns there).
/// A knot that leaves a spline only continuous by its multiplicity, where the spline's
/// derivative is the same either side (`eval`, the curve's own evaluator), is once
/// differentiable and gets a pcurve with a tangent there: two pieces fitted apart would meet at
/// an angle, sharp in (u, v) where the surface's parameters run at very different speeds (a
/// thread's flank, hundreds of mm per unit along the helix and one across it).
fn breaks(curve: &Curve, eval: Option<&CurveEval>, t0: f64, t1: f64) -> Vec<(f64, usize)> {
let spline = |knots: &[f64], count: usize| -> Vec<(f64, usize)> {
let degree = knots.len().saturating_sub(count + 1);
let mut out: Vec<(f64, usize)> = Vec::new();
for &k in knots.iter().filter(|&&k| k > t0 && k < t1) {
match out.last_mut() {
Some(last) if last.0 == k => last.1 += 1,
_ => out.push((k, 1)),
}
}
out.into_iter()
.map(|(k, m)| {
let mut smooth = degree.saturating_sub(m);
if m == degree && eval.is_some_and(|e| differentiable_at(e, k)) {
smooth = 1;
}
(k, (3 - smooth.min(2)).clamp(1, 3))
})
.collect()
};
let corners = |n: usize| (1..n.saturating_sub(1)).map(|i| (i as f64, 3)).filter(|k| k.0 > t0 && k.0 < t1).collect();
match curve {
Curve::Line { .. } | Curve::Circle { .. } | Curve::Ellipse { .. } => Vec::new(),
Curve::Polyline { points } => corners(points.len()),
Curve::Spline { knots, points } => spline(knots, points.len()),
Curve::RationalSpline { knots, homogeneous } => spline(knots, homogeneous.len()),
Curve::SurfaceIntersection { leader, .. } => match &**leader {
Curve::Polyline { points } => corners(points.len()),
other => breaks(other, None, t0, t1),
},
}
}
/// Relative difference of a spline's derivatives either side of a knot within which it is
/// differentiable there (rounding of derivatives made equal).
const SAME_DERIVATIVE: f64 = 1e-9;
/// Whether the spline `eval` has one derivative at its knot `k` (the same from either side).
fn differentiable_at(eval: &CurveEval, k: f64) -> bool {
let CurveEval::Spline(c) = eval else { return false };
let Some((l, r)) = c.sides(k) else { return false };
let size = l.len().max(r.len());
size > 0.0 && (l - r).len() <= SAME_DERIVATIVE * size
}
/// Signed area in (u, v) of the loop with its uses moved by `shift` along `d`.
fn area(uses: &[(Curve2, bool)], shift: &[f64], d: usize) -> f64 {
let mut ring: Vec<[f64; 2]> = Vec::new();
@@ -643,127 +498,6 @@ fn area(uses: &[(Curve2, bool)], shift: &[f64], d: usize) -> f64 {
(0..n).filter_map(|i| Some((ring.get(i)?, ring.get((i + 1) % n)?))).map(|(p, q)| p[0] * q[1] - q[0] * p[1]).sum::<f64>() * 0.5
}
/// The samples' parameters, where a pole leaves one undefined taken from the nearest sample
/// where it is defined (so the pcurve runs straight into the pole, never across to another
/// angle).
fn fill_poles(samples: &[Sample]) -> Vec<(f64, [f64; 2])> {
let mut out: Vec<(f64, [f64; 2])> = samples.iter().map(|s| (s.t, s.uv)).collect();
for d in 0..2 {
// Where the surface collapses along `d`, the parameter `d` says nothing.
let free = |s: &Sample| s.pole.get(d).copied().unwrap_or(false);
if samples.iter().all(free) {
continue;
}
for (i, s) in samples.iter().enumerate() {
if !free(s) {
continue;
}
let left = samples.get(..i).and_then(|l| l.iter().rposition(|x| !free(x)));
let right = samples.get(i + 1..).and_then(|r| r.iter().position(|x| !free(x))).map(|k| k + i + 1);
let dt = |k: usize| samples.get(k).map_or(f64::INFINITY, |x| (x.t - s.t).abs());
let pick = match (left, right) {
(Some(l), Some(r)) => Some(if dt(l) <= dt(r) { l } else { r }),
(a, b) => a.or(b),
};
if let (Some(src), Some(slot)) = (pick.and_then(|k| samples.get(k)), out.get_mut(i))
&& let (Some(x), Some(v)) = (slot.1.get_mut(d), src.uv.get(d))
{
*x = *v;
}
}
}
out
}
/// Cubic control points on `knots` through the first and last data points, least squares
/// through the rest (banded normal equations, Cholesky).
fn least_squares(knots: &[f64], data: &[(f64, [f64; 2])]) -> Option<Vec<[f64; 2]>> {
const P: usize = 3;
let n = knots.len().checked_sub(P + 1)?;
let (first, last) = (data.first()?.1, data.last()?.1);
let m = n.checked_sub(2)?;
// Normal equations, upper band: ata[j][k] = A(j, j + k) over the free points 1..n−1.
let mut ata = vec![[0.0f64; P + 1]; m];
let mut atb = vec![[0.0f64; 2]; m];
for &(t, q) in data {
let span = crate::bspline::find_span(knots, P, n, t)?;
let nb = crate::tess::basis_ders(knots, span, P, t)?;
let nb = nb.first()?.get(..=P)?;
let mut r = q;
for (j, w) in nb.iter().enumerate() {
let i = span + j - P;
let fixed = if i == 0 {
Some(first)
} else if i + 1 == n {
Some(last)
} else {
None
};
if let Some(f) = fixed {
r = [r[0] - w * f[0], r[1] - w * f[1]];
}
}
for (j, wj) in nb.iter().enumerate() {
let Some(row) = (span + j - P).checked_sub(1).filter(|&r| r < m) else { continue };
if let Some(b) = atb.get_mut(row) {
b[0] += wj * r[0];
b[1] += wj * r[1];
}
for (k, wk) in nb.iter().enumerate().skip(j) {
let Some(col) = (span + k - P).checked_sub(1).filter(|&c| c < m) else { continue };
if let Some(x) = ata.get_mut(row).and_then(|a| a.get_mut(col - row)) {
*x += wj * wk;
}
}
}
}
// A = Rᵀ R with R upper of the same band: r[j][k] = R(j, j + k).
let at = |r: &[[f64; P + 1]], j: usize, k: usize| r.get(j).and_then(|x| x.get(k)).copied().unwrap_or(0.0);
let scale = ata.iter().map(|a| a[0]).fold(0.0, f64::max);
let mut r = vec![[0.0f64; P + 1]; m];
for j in 0..m {
let mut diag = at(&ata, j, 0);
for l in 1..=P.min(j) {
diag -= at(&r, j - l, l).powi(2);
}
if !(diag > 1e-14 * scale) {
return None;
}
let diag = diag.sqrt();
let mut row = [diag, 0.0, 0.0, 0.0];
for (k, x) in row.iter_mut().enumerate().skip(1).take_while(|(k, _)| j + k < m) {
let mut v = at(&ata, j, k);
for l in 1..=(P - k).min(j) {
v -= at(&r, j - l, l) * at(&r, j - l, l + k);
}
*x = v / diag;
}
*r.get_mut(j)? = row;
}
let mut y = vec![[0.0f64; 2]; m];
for j in 0..m {
let mut v = *atb.get(j)?;
for l in 1..=P.min(j) {
let (c, yl) = (at(&r, j - l, l), *y.get(j - l)?);
v = [v[0] - c * yl[0], v[1] - c * yl[1]];
}
let d = at(&r, j, 0);
*y.get_mut(j)? = [v[0] / d, v[1] / d];
}
let mut x = vec![[0.0f64; 2]; m];
for j in (0..m).rev() {
let mut v = *y.get(j)?;
for k in (1..=P).take_while(|k| j + k < m) {
let (c, xk) = (at(&r, j, k), *x.get(j + k)?);
v = [v[0] - c * xk[0], v[1] - c * xk[1]];
}
let d = at(&r, j, 0);
*x.get_mut(j)? = [v[0] / d, v[1] / d];
}
let out: Vec<[f64; 2]> = std::iter::once(first).chain(x).chain(std::iter::once(last)).collect();
out.iter().flatten().all(|v| v.is_finite()).then_some(out)
}
/// The control points of `c` in (u, v).
fn points(c: &Curve2) -> Vec<[f64; 2]> {
match c {
crates/step/src/topology.rs
+21 -44
View File
@@ -694,17 +694,14 @@ impl<'b, 'a> Builder<'b, 'a> {
let mut face_forms: Vec<Option<BSurface>> = Vec::with_capacity(self.faces.len());
let mut domains = Vec::with_capacity(self.faces.len());
let mut turned = vec![false; self.faces.len()];
let (mut failed, mut first_failure, mut loosened) = (0usize, None, 0usize);
let (mut failed, mut first_failure) = (0usize, None);
for fi in 0..self.faces.len() {
match self.face_pcurves(fi) {
Ok((pc, domain, flipped, tol, form)) => {
Ok((pc, domain, flipped, form)) => {
face_forms.push(form);
if flipped && let Some(t) = turned.get_mut(fi) {
*t = true;
}
if tol > self.tol {
loosened += 1;
}
pcs.push(pc);
domains.push(domain);
}
@@ -724,9 +721,6 @@ impl<'b, 'a> Builder<'b, 'a> {
if flipped > 0 {
notes.push(format!("{flipped} faces turned over: their loops disagree with their same_sense"));
}
if loosened > 0 {
notes.push(format!("{loosened} faces have their pcurves fitted looser than {} mm", self.tol));
}
let shells: Vec<Vec<(u32, bool)>> = shells
.into_iter()
.map(|s| s.into_iter().map(|(f, fwd)| (f, if turned.get(f as usize).copied().unwrap_or(false) { !fwd } else { fwd })).collect())
@@ -953,11 +947,11 @@ impl<'b, 'a> Builder<'b, 'a> {
}
/// The pcurves of a face's edge loops (in its edge-loop order), its domain, whether it was
/// turned over, and the tolerance they were fitted at. Picks the outer loop when the file
/// turned over, and its NURBS form. Picks the outer loop when the file
/// marks none (or marks one that another loop encloses), and turns the face over when its
/// outer loop runs clockwise about the natural normal on a surface where that cannot be a
/// face's outer boundary.
fn face_pcurves(&mut self, fi: usize) -> R<(Vec<Vec<Curve2>>, [[f64; 2]; 2], bool, f64, Option<BSurface>)> {
fn face_pcurves(&mut self, fi: usize) -> R<(Vec<Vec<Curve2>>, [[f64; 2]; 2], bool, Option<BSurface>)> {
let face = self.faces.get(fi).ok_or("face reference")?;
let label = face.surface.clone();
let edge_loops: Vec<Vec<(u32, bool)>> = face
@@ -970,52 +964,35 @@ impl<'b, 'a> Builder<'b, 'a> {
.collect();
// Which loop index in `face.loops` each edge loop is.
let positions: Vec<usize> = face.loops.iter().enumerate().filter(|(_, l)| matches!(l, LoopB::Edges(_))).map(|(i, _)| i).collect();
let mut tol = self.tol;
let tol = self.tol;
// The face's pcurves on `on`, each edge on its NURBS form when `formed`.
let pcurves_on = |on: &BSurface, formed: bool, tol: &mut f64| -> R<Vec<Vec<Curve2>>> {
loop {
let mut uses = Vec::with_capacity(edge_loops.len());
for l in &edge_loops {
let mut u = Vec::with_capacity(l.len());
for &(e, forward) in l {
let edge = self.edges.get(e as usize).ok_or("edge reference")?;
let (curve, interval) = match &edge.form {
Some((c, i)) if formed => (c, *i),
_ => (&edge.curve, edge.interval),
};
u.push(EdgeUse { edge: e, curve, interval, forward, tolerance: 0.0, given: None });
}
uses.push(u);
}
match pcurve::pcurves(on, &uses, *tol) {
Ok(p) => return Ok(p),
Err(pcurve::Error::Failed(m)) if !m.contains("not supported") && *tol * 10.0 <= MAX_FIT_TOLERANCE * (1.0 + 1e-9) => {
*tol *= 10.0;
}
Err(e) => return Err(e.to_string()),
let pcurves_on = |on: &BSurface, formed: bool| -> R<Vec<Vec<Curve2>>> {
let mut uses = Vec::with_capacity(edge_loops.len());
for l in &edge_loops {
let mut u = Vec::with_capacity(l.len());
for &(e, forward) in l {
let edge = self.edges.get(e as usize).ok_or("edge reference")?;
let (curve, interval) = match &edge.form {
Some((c, i)) if formed => (c, *i),
_ => (&edge.curve, edge.interval),
};
u.push(EdgeUse { edge: e, curve, interval, forward, tolerance: 0.0, given: None });
}
uses.push(u);
}
pcurve::pcurves(on, &uses, tol).map_err(|e| e.to_string())
};
// A cylinder, cone, sphere or torus on its NURBS form, over the extent its label's pcurves
// give; then every pcurve on the forms.
let form = match brep::form::frame_of(&label) {
Some(_) => {
let extent = face_domain(&label, &pcurves_on(&label, false, &mut tol)?)?;
let extent = face_domain(&label, &pcurves_on(&label, false)?)?;
Some(brep::form::surface(&label, extent).ok_or_else(|| unsupported("a face of more than one turn"))?.0)
}
None => None,
};
let surface = form.clone().unwrap_or(label);
let pcs = pcurves_on(&surface, true, &mut tol)?;
if tol > self.tol {
for l in &edge_loops {
for &(e, _) in l {
if let Some(edge) = self.edges.get_mut(e as usize) {
edge.tolerance = edge.tolerance.max(tol);
}
}
}
}
let pcs = pcurves_on(&surface, true)?;
// Signed areas of the loops that close in (u, v).
let periods = natural_periods(&surface);
// Each use starts near where the one before it ends, along a bounded direction that
@@ -1077,7 +1054,7 @@ impl<'b, 'a> Builder<'b, 'a> {
}
let pcs: Vec<Vec<Curve2>> = if flip { pcs.into_iter().map(|l| l.into_iter().rev().collect()).collect() } else { pcs };
let domain = face_domain(&surface, &pcs)?;
Ok((pcs, domain, flip, tol, form))
Ok((pcs, domain, flip, form))
}
}