brep: a seam from pole to pole gets pcurves a period apart

A whole sphere bounded only by its seam (SMLib's sphere: one meridian
used twice) came out with both uses on the same u line, a loop that
encloses nothing. Two causes in the per-loop closure:

- The pole-junction search took any positive signed area as the right
  turn. With the two uses on one line the area is rounding, and in most
  frames that rounding is a tiny positive number, so "no shift" won over
  the period shift. A turn now needs more than 1e-9 of the loop's
  period box.
- Junctions counted as poles only where the surface collapses exactly.
  Ends within the tolerance of a pole (as a kernel's records have them)
  were closed to the nearest period instead. A junction is now free
  where moving its periodic parameter by half a period stays within the
  tolerance, and fit samples within a quarter of it take their free
  parameter from a neighbour, so the pcurve runs straight into the pole.

Tests: the sphere bounded by its seam in four frames, three meridian
angles, circle and NURBS meridians, both loop orders, ends at or a tenth
of the tolerance off the poles; and a torus bounded by its two seams.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01QWhMwdqJYQKQq8RCLSL1LT
This commit is contained in:
jyviko
2026-10-21 01:00:51 -04:00
parent 1618329763
commit 25d2056a2e
2 changed files with 126 additions and 11 deletions
+39 -11
View File
@@ -101,7 +101,7 @@ pub fn pcurves(surface_record: &Surface, loops: &[Vec<EdgeUse<'_>>], tolerance:
};
fitted.push((pc, u.forward));
}
on.close(&mut fitted)?;
on.close(&mut fitted, tolerance)?;
out.push(fitted.into_iter().map(|(c, _)| c).collect());
}
Ok(out)
@@ -254,10 +254,29 @@ impl On {
Some(q)
}
fn sample(&self, curve: &CurveEval, t: f64, near: Option<[f64; 2]>, accept: f64) -> Result<Sample> {
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.surface.pole(uv), mid: None })
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.
fn pole(&self, uv: [f64; 2], reach: f64) -> [bool; 2] {
let mut out = self.surface.pole(uv);
let Some((_, su, sv)) = self.surface.derivs(uv[0], uv[1]) else { return out };
for (d, speed) in [su.len(), sv.len()].into_iter().enumerate() {
let span = match (self.dom.period.get(d).copied().flatten(), self.dom.range.get(d).copied().flatten()) {
(Some(period), _) => 0.5 * period,
(None, Some((a, b))) => (b - a).abs(),
(None, None) => continue,
};
if let Some(flag) = out.get_mut(d) {
*flag |= speed * span <= reach;
}
}
out
}
/// The error for a curve point `p` (at `t`) farther than the fit allows from the surface.
@@ -272,6 +291,9 @@ impl On {
/// 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, t0, t1);
@@ -289,7 +311,7 @@ impl On {
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)?;
let s = self.sample(curve, t, near, accept, reach)?;
near = Some(s.uv);
samples.push(s);
}
@@ -340,7 +362,7 @@ impl On {
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.surface.pole(uv), mid: None });
dense.push(Sample { t, p: m, uv, pole: self.pole(uv, reach), mid: None });
}
_ => dense.push(s),
}
@@ -352,10 +374,11 @@ impl On {
/// 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 the periodic parameter is free:
/// there the offsets that turn the loop counter-clockwise in (u, v) (a face boundary seen
/// from the natural normal) are taken, else the nearest.
fn close(&self, uses: &mut [(Curve2, bool)]) -> Result<()> {
/// loop's lowest point allows. Where uses meet at a pole (within `tol`) the periodic
/// parameter is free: there the offsets that turn the loop counter-clockwise in (u, v) (a
/// face boundary seen from the natural normal) by more than rounding are taken, else the
/// nearest. A seam used twice from pole to pole thus encloses its face, never nothing.
fn close(&self, uses: &mut [(Curve2, bool)], tol: f64) -> Result<()> {
let n = uses.len();
let mut ends = Vec::with_capacity(n);
for (c, fwd) in uses.iter() {
@@ -367,7 +390,7 @@ impl On {
let Some(period) = self.dom.period.get(d).copied().flatten().filter(|p| *p > 0.0) else { continue };
// Start after a pole junction if there is one: from there continuity fixes the
// offsets up to the next pole.
let at_pole: Vec<bool> = ends.iter().map(|(s, _)| self.surface.pole(*s).get(d).copied().unwrap_or(false)).collect();
let at_pole: Vec<bool> = ends.iter().map(|(s, _)| self.pole(*s, tol).get(d).copied().unwrap_or(false)).collect();
let first = at_pole.iter().position(|&p| p).unwrap_or(0);
let order: Vec<usize> = (0..n).map(|k| (first + k) % n).collect();
let mut shift = vec![0.0; n];
@@ -384,6 +407,11 @@ impl On {
// After each later pole one period less, none or more: the right turn, then the
// fewest moves.
if !free.is_empty() && free.len() <= 6 {
// Less than this is rounding: a loop that encloses nothing (a seam's two uses
// on one line) has no turn, whatever the sign of its area.
let other = uses.iter().flat_map(|(c, _)| points(c).into_iter().filter_map(|q| q.get(1 - d).copied()));
let (lo, hi) = other.fold((f64::INFINITY, f64::NEG_INFINITY), |(a, b), x| (a.min(x), b.max(x)));
let flat = if hi >= lo { 1e-9 * period * (hi - lo) } else { 0.0 };
let mut best: Option<(bool, usize, Vec<f64>)> = None;
for combo in 0..3usize.pow(free.len() as u32) {
let (mut code, mut moved, mut add) = (combo, 0, 0.0);
@@ -399,7 +427,7 @@ impl On {
*slot += add;
}
}
let right = area(uses, &extra, d) > 0.0;
let right = area(uses, &extra, d) > flat;
if best.as_ref().is_none_or(|(r, m, _)| (right && !r) || (right == *r && moved < *m)) {
best = Some((right, moved, extra));
}
+87
View File
@@ -234,6 +234,93 @@ fn meridians_run_straight_into_the_poles() {
}
}
/// A whole sphere bounded only by its seam (one meridian from pole to pole, used twice), in
/// frames where rounding leaves the loop with the uses on one line a tiny positive area: the
/// two uses come out a period apart and the loop turns counter-clockwise around the whole
/// (u, v) box, in either loop order, for a circle or a NURBS meridian, with its ends at the
/// poles or within a tenth of the tolerance of them.
#[test]
fn a_seam_from_pole_to_pole_encloses_the_whole_sphere() {
let w = std::f64::consts::FRAC_1_SQRT_2;
let tol = 1e-7;
let frames = [
about_z([0.; 3], false),
about_z([1., -2., 3.], true),
frame([1., 2., 3.], [0., 0.6, 0.8], [0., -0.8, 0.6]),
frame([-4., 0.5, 2.], [0.36, 0.48, 0.8], [-0.8, 0.6, 0.]),
];
for sf in frames {
let sphere = Surface::Sphere { frame: sf.clone(), radius: 6. };
let ax = crate::brep::eval::Axes::new(&sf);
for turn in [0f64, 2.1, -0.7] {
let x = ax.x * turn.cos() + ax.y * turn.sin();
let p = |v: Vec3| [v.x, v.y, v.z];
let (o, r, z) = (ax.o, x * 6., ax.z * 6.);
let circle = Curve::Circle { frame: Frame { origin: p(o), x_axis: p(x), y_axis: p(ax.z) }, radius: 6. };
for off in [0., 0.1 * tol] {
// South pole to north pole through the equator, as two rational quarter arcs,
// the ends `off` to the side of the poles.
let h = |v: Vec3, wt: f64| [v.x * wt, v.y * wt, v.z * wt, wt];
let side = ax.y * off;
let nurbs = Curve::RationalSpline {
knots: vec![0., 0., 0., 0.5, 0.5, 1., 1., 1.],
homogeneous: vec![h(o - z + side, 1.), h(o - z + r, w), h(o + r, 1.), h(o + z + r, w), h(o + z - side, 1.)],
};
for (seam, iv) in [(&circle, [-FRAC_PI_2, FRAC_PI_2]), (&nurbs, [0., 1.])] {
for order in [[true, false], [false, true]] {
let l: Vec<EdgeUse> = order.iter().map(|&forward| EdgeUse { edge: 0, curve: seam, interval: iv, forward }).collect();
let pcs = &pcurves(&sphere, std::slice::from_ref(&l), tol).unwrap()[0];
let what = format!("{sf:?} {turn} {off} {seam:?} {order:?}");
for pc in pcs {
assert!(lift_error(&sphere, seam, iv, pc) <= tol, "{what}");
}
let mid = 0.5 * (iv[0] + iv[1]);
let (a, b) = (pcs[0].eval(mid).unwrap(), pcs[1].eval(mid).unwrap());
assert!(((a[0] - b[0]).abs() - TAU).abs() < 1e-6 && (a[1] - b[1]).abs() < 1e-12, "a period apart: {a:?} {b:?} {what}");
let (ends, area) = walk(&l, pcs);
for k in 0..2 {
assert!((ends[k].1[1] - ends[(k + 1) % 2].0[1]).abs() < 1e-6, "closed in v: {ends:?} {what}");
}
// The box less the corners `walk` cuts (it joins each use's 16th
// sample to the next use's start): about 15/16 of it.
assert!((area / (TAU * PI * 15. / 16.) - 1.).abs() < 1e-2, "counter-clockwise around the box: {area} {what}");
}
}
}
}
}
}
/// A whole torus bounded only by its two seams (a circle about the axis and a circle around
/// the tube, each used twice): no pole, continuity alone puts each seam's uses a period apart
/// and the loop around the whole (u, v) box.
#[test]
fn a_torus_bounded_by_its_seams_spans_both_periods() {
let frames = [about_z([0.; 3], false), frame([1., 2., 3.], [0., 0.6, 0.8], [0., -0.8, 0.6])];
for sf in frames {
let torus = Surface::Torus { frame: sf.clone(), major_radius: 10., minor_radius: 2. };
let ax = crate::brep::eval::Axes::new(&sf);
let p = |v: Vec3| [v.x, v.y, v.z];
// The outer equator (v = 0, u from 0) and the tube circle at u = 0 (v from 0, upwards).
let around = Curve::Circle { frame: sf.clone(), radius: 12. };
let tube = Curve::Circle { frame: Frame { origin: p(ax.o + ax.x * 10.), x_axis: p(ax.x), y_axis: p(ax.z) }, radius: 2. };
let full = [0., TAU];
let l = vec![
EdgeUse { edge: 0, curve: &around, interval: full, forward: true },
EdgeUse { edge: 1, curve: &tube, interval: full, forward: true },
EdgeUse { edge: 0, curve: &around, interval: full, forward: false },
EdgeUse { edge: 1, curve: &tube, interval: full, forward: false },
];
let pcs = &pcurves(&torus, std::slice::from_ref(&l), 1e-9).unwrap()[0];
let (ends, area) = walk(&l, pcs);
for k in 0..4 {
let (end, start) = (ends[k].1, ends[(k + 1) % 4].0);
assert!((end[0] - start[0]).abs() < 1e-12 && (end[1] - start[1]).abs() < 1e-12, "{k}: {end:?} {start:?}");
}
assert!((area / (TAU * TAU) - 1.).abs() < 1e-12, "counter-clockwise around the box: {area}");
}
}
/// An intersection record evaluates to the point of both surfaces in the leader's normal
/// plane, and its pcurves on either surface fit.
#[test]