From 7b79e7cb617258765192e77bb2e45fe1fa882662 Mon Sep 17 00:00:00 2001 From: Akhil Velagapudi Date: Sun, 11 Oct 2026 06:32:43 +0000 Subject: [PATCH] Bound spline patches in a frame on each parameter direction (ENG-3833) Projecting a point onto a spline surface refines Bezier patches until none can be closer than the incumbent. One lower bound measures each patch in the frame (normal, S_u, normal x S_u). When S_v is oblique to S_u, a long, thin patch extends along S_u by its length times the skew, not by its width, and the bound misses the offset of a point beside it. RFPC-SMA16-FN-A has a spline face with a knot span 2e-4 mm wide next to spans 2.6 mm long, with S_v 3.8% oblique. Its patches across the thin span could not be told apart, so each split of the span doubled the live patches: some boundary points examined 50,000 patches, and two faces took 1.6 s and 1 s on one core. The patch is now also bounded in the frame (normal, S_v, normal x S_v), and the larger bound is kept. Both are valid bounds, so the search returns the same minimum to within roundoff. RFPC-SMA16-FN-A now takes 0.3 s instead of 1.8 s. --- nurbs/src/sampled_surface.rs | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/nurbs/src/sampled_surface.rs b/nurbs/src/sampled_surface.rs index c35da05..d981246 100644 --- a/nurbs/src/sampled_surface.rs +++ b/nurbs/src/sampled_surface.rs @@ -563,16 +563,25 @@ where let d = self.surf.derivs_in_span::<2>(mid, spans, p); let normal = d[1][0].cross(&d[0][1]); let normal = if normal.norm_squared() > 0. { normal.normalize() } else { DVec3::zeros() }; - let tangent = if d[1][0].norm_squared() > 0. { d[1][0].normalize() } else { DVec3::zeros() }; - let bitangent = normal.cross(&tangent); - let mut slab = [DVec3::repeat(f64::INFINITY), DVec3::repeat(f64::NEG_INFINITY)]; + // Parameter directions need not be orthogonal. A frame built on + // one of them measures the other's skew as extent across it: a + // long, thin patch then spans its own width many times over and + // hides the offset of every point beside it. Bound in a frame on + // each direction, so the long one carries the skew. + let frames = [d[1][0], d[0][1]].map(|tangent| { + let tangent = if tangent.norm_squared() > 0. { tangent.normalize() } else { DVec3::zeros() }; + [normal, tangent, normal.cross(&tangent)] + }); + let mut slabs = [[DVec3::repeat(f64::INFINITY), DVec3::repeat(f64::NEG_INFINITY)]; 2]; let mut bounds = [DVec3::repeat(f64::INFINITY), DVec3::repeat(f64::NEG_INFINITY)]; let direction = residual.unwrap_or_else(DVec3::zeros).try_normalize(0.).unwrap_or_else(DVec3::zeros); let mut support = f64::INFINITY; for q in patches[patch].iter().copied().map(crate::nd_curve::cartesian) { bounds[0] = bounds[0].inf(&q); bounds[1] = bounds[1].sup(&q); - let offset = DVec3::new((q-p).dot(&normal), (q-p).dot(&tangent), (q-p).dot(&bitangent)); - slab[0] = slab[0].inf(&offset); slab[1] = slab[1].sup(&offset); + for (slab, frame) in slabs.iter_mut().zip(&frames) { + let offset = DVec3::from_fn(|i, _| (q-p).dot(&frame[i])); + slab[0] = slab[0].inf(&offset); slab[1] = slab[1].sup(&offset); + } support = support.min((q-p).dot(&direction)); } // Any supporting plane bounds the hull. The incumbent residual @@ -580,17 +589,17 @@ where // lose, especially at constrained endpoint minima. // An orthonormal control box also bounds distance beyond a finite // patch edge; the normal slab alone only bounds its infinite plane. - let slab_roundoff = DVec3::new(normal.abs().dot(&roundoff), - tangent.abs().dot(&roundoff), bitangent.abs().dot(&roundoff)); + let slab_roundoff = frames.map(|frame| DVec3::from_fn(|i, _| frame[i].abs().dot(&roundoff))); // Stop when distance comparisons overlap at input resolution, // not when an optimistic bound beats an optimistic incumbent. // Fixed arithmetic uncertainty cannot shrink with subdivision. // This is a numerical tie, not certified exact domination. let box_gap = (bounds[0]-p).sup(&(p-bounds[1])); - let slab_gap = slab[0].sup(&(-slab[1])); + let slab_gaps = slabs.map(|slab| slab[0].sup(&(-slab[1]))); let support_roundoff = direction.abs().dot(&roundoff); - let bound = |sign: f64| (box_gap + sign * roundoff).sup(&DVec3::zeros()).norm() - .max((slab_gap + sign * slab_roundoff).sup(&DVec3::zeros()).norm()) + let bound = |sign: f64| (0..2).map(|i| (slab_gaps[i] + sign * slab_roundoff[i]) + .sup(&DVec3::zeros()).norm()) + .fold((box_gap + sign * roundoff).sup(&DVec3::zeros()).norm(), f64::max) .max(support + sign * support_roundoff); // Definitely farther patches need no candidate evaluations. if residual.is_some_and(|r| bound(-1.) >= (r.abs() + roundoff).norm()) {