brep: the projection runs Newton from up to 32 starts, nearest first, until one lands

On a thread's root band (1.4 mm tall at a pitch of 1.5) a point on the band's lower rim lies
0.1 mm from the rim of the turn below. The search's starts are dense across the band and
sparse along it, so its four nearest starts all lay on the turn below, Gauss-Newton settled
there, and the pcurve fit refused the edge ("the curve at 0 lies 0.1008 mm off the surface",
ladder/06-threads/007). The search now goes on through the nearest 32 starts until one
converges within the accepted distance, and stops at the first that does.

The test strips now put their knots at other angles on each turn (7.3 control points a
turn), as SMLib's screw sweeps do.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017itGh5fpMEousnZ91cYcVt
This commit is contained in:
jyviko
2026-10-21 15:23:31 -04:00
parent e1c631dd33
commit 829d161f6b
2 changed files with 41 additions and 8 deletions
+32 -7
View File
@@ -486,22 +486,30 @@ fn intersection_curves_lie_on_both_surfaces() {
}
}
/// A thread flank as a screw sweep leaves it: a spline strip winding twenty turns of pitch 1 about
/// A thread flank as a screw sweep leaves it: a spline strip winding twenty-two turns of pitch 1 about
/// z, rising from radius 3 to radius 4 by half the pitch, so one turn lies 0.89 mm from the next,
/// and its v = ½ line (its rows' midpoints on its u knots, so it lies on the strip).
fn thread_flank() -> (Surface, Curve) {
const TURNS: usize = 20;
const PER_TURN: usize = 8;
let n = TURNS * PER_TURN + 1;
helical_strip([3., 4.], 0.5, 0.5)
}
/// A spline strip winding twenty-two turns of pitch 1 about z from radius `r[0]` (v = 0) to `r[1]`
/// (v = 1), rising by `rise` across, and its line at v = `at` (its rows' points there on its u
/// knots, so it lies on the strip).
fn helical_strip(r: [f64; 2], rise: f64, at: f64) -> (Surface, Curve) {
// Control points 7.3 to a turn (161 over 22 turns), so the knots fall at other angles on
// each turn.
const PER_TURN: f64 = 7.3;
let n = 161;
let helix = |r: f64, rise: f64| -> Vec<[f64; 3]> {
(0..n)
.map(|i| {
let a = TAU * i as f64 / PER_TURN as f64;
let a = TAU * i as f64 / PER_TURN;
[r * a.cos(), r * a.sin(), a / TAU + rise]
})
.collect()
};
let (inner, outer) = (helix(3., 0.), helix(4., 0.5));
let (inner, outer) = (helix(r[0], 0.), helix(r[1], rise));
let mut u_knots = vec![0.; 4];
u_knots.extend((1..n - 3).map(|k| k as f64 / (n - 3) as f64));
u_knots.extend([1.; 4]);
@@ -510,7 +518,8 @@ fn thread_flank() -> (Surface, Curve) {
v_knots: vec![0., 0., 1., 1.],
points: inner.iter().zip(&outer).map(|(a, b)| vec![*a, *b]).collect(),
};
let middle = Curve::Spline { knots: u_knots, points: inner.iter().zip(&outer).map(|(a, b)| [0, 1, 2].map(|k| 0.5 * (a[k] + b[k]))).collect() };
let middle =
Curve::Spline { knots: u_knots, points: inner.iter().zip(&outer).map(|(a, b)| [0, 1, 2].map(|k| (1. - at) * a[k] + at * b[k])).collect() };
(surface, middle)
}
@@ -528,6 +537,22 @@ fn edges_on_a_spline_winding_several_turns_stay_on_their_turn() {
}
}
/// A thread's root band, 0.93 of the pitch tall, its turns 0.07 mm apart: an edge along its lower
/// rim (as where the root meets a shoulder) gets a pcurve on its own turn, though the rim of the
/// turn below passes nearer some of the search's starts than its own turn does.
#[test]
fn an_edge_along_the_rim_of_a_band_winding_close_to_itself_stays_on_its_turn() {
let (surface, rim) = helical_strip([5., 5.], 0.93, 0.);
// Short edges starting every 1/64 of a turn round one turn.
for k in 0..64 {
let t0 = 0.5 + f64::from(k) / (64. * 22.);
let interval = [t0, t0 + 1. / 64.];
let pc = &pcurves(&surface, &one(&rim, interval), 2e-4).unwrap()[0][0];
let e = lift_error(&surface, &rim, interval, pc);
assert!(e <= 2e-4, "{interval:?}: {e}");
}
}
/// A curve that is smooth through a knot of full multiplicity (one smooth curve with the knot
/// inserted three times, as SMLib's intersection curves come) gets a pcurve with one tangent
/// there. Two pieces fitted apart meet at an angle in (u, v) (a tenth of a degree here, over the
+9 -1
View File
@@ -210,12 +210,18 @@ pub fn project(s: &EvalSurface, dom: &Domain, p: Vec3, hint: Option<[f64; 2]>, a
}
}
starts.sort_by(|a, b| a.0.total_cmp(&b.0));
for (_, start) in starts.iter().take(4) {
// The nearest starts first, until one converges within `accept`: where turns of a surface
// pass close to each other (a thread's root band a tenth of a millimetre from the next
// turn's), the nearest few can all lie on the wrong turn.
for (_, start) in starts.iter().take(NEWTON_STARTS) {
if let Some((uv, d)) = newton(s, dom, p, *start)
&& best.is_none_or(|(_, bd)| d < bd)
{
best = Some((uv, d));
}
if accept.is_finite() && best.is_some_and(|(_, d)| d <= accept) {
break;
}
}
best.filter(|(_, d)| *d <= accept).map(|(uv, _)| uv)
}
@@ -226,6 +232,8 @@ const EVEN_STARTS: usize = 16;
const SPAN_STARTS: usize = 2;
/// Most points of the surface the search measures.
const MAX_STARTS: usize = 1 << 14;
/// Most of them it runs Gauss–Newton from.
const NEWTON_STARTS: usize = 32;
/// The parameters along direction `d` (over `range`) that [`project`] starts from when it has no
/// hint: even steps across the range, and where the surface has its own knots along `d` (a