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

brep: a surface point's parameters move to geom

Browse Source 
The parameter box, the projection of a point onto an evaluator surface (closed form on a
plane, damped Gauss-Newton elsewhere), unwrapping and the pole test work on geom's
EvalSurface alone. They move from the kernel's tessellator to solidcraft_geom::tess::uv so a
pcurve routine on plain B-rep records can use them without truck; the kernel re-exports them
and its tessellation is unchanged.

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-20 23:49:36 -04:00
parent a8c8f88323
commit 1d49160081
3 changed files with 190 additions and 188 deletions
Download Patch File Download Diff File Expand all files Collapse all files
crates/geom/src/tess.rs
+4 -1
View File
@@ -9,11 +9,14 @@
//! [`EvalSurface`] forms defined here, so the two cannot drift apart.
//!
//! References: Piegl & Tiller, *The NURBS Book* (2nd ed.), A2.1 (span), A2.3 (basis
//! functions and derivatives), §4.5 (derivatives of a rational surface).
//! functions and derivatives), §4.5 (derivatives of a rational surface). [`uv`] finds the
//! parameters of a point on an evaluator surface.
use crate::bspline::MAX_SPLINE_DEGREE;
use crate::{Mesh, Vec3};
pub mod uv;
/// Most knots in one evaluable knot vector (bounds the shader's span search).
pub const MAX_EVAL_KNOTS: usize = 1 << 12;
/// Most control points in one evaluable surface.
crates/geom/src/tess/uv.rs
+184
View File
@@ -0,0 +1,184 @@
//! Parameters of points on an [`EvalSurface`], on the evaluator's own surface form: a
//! surface's parameter box and periods, the parameters of a point on it (closed form on a plane,
//! damped Gauss–Newton elsewhere), unwrapping across a period, and its poles.
use super::EvalSurface;
use crate::Vec3;
/// The parameter box of a surface and its periods.
#[derive(Clone, Copy, Debug)]
pub struct Domain {
/// Bounded ranges (`None`: unbounded, as a plane's).
pub range: [Option<(f64, f64)>; 2],
pub period: [Option<f64>; 2],
}
/// The parameter box and periods of `s`.
pub fn domain(s: &EvalSurface) -> Domain {
let tau = std::f64::consts::TAU;
match s {
EvalSurface::Plane { .. } => Domain { range: [None, None], period: [None, None] },
EvalSurface::Nurbs { degree, knots, count, .. } => {
let r = |d: usize| Some((*knots[d].get(degree[d])?, *knots[d].get(count[d])?));
Domain { range: [r(0), r(1)], period: [None, None] }
}
EvalSurface::Revolved { profile, swap, .. } => {
let r = (|| Some((*profile.knots.get(profile.degree)?, *profile.knots.get(profile.ctrl.len())?)))();
let (a, b) = ((r, None), (Some((0.0, tau)), Some(tau)));
let (u, v) = if *swap { (b, a) } else { (a, b) };
Domain { range: [u.0, v.0], period: [u.1, v.1] }
}
}
}
impl Domain {
pub fn clamp(&self, uv: [f64; 2]) -> [f64; 2] {
let mut out = uv;
for (d, x) in out.iter_mut().enumerate() {
if let (Some((a, b)), None) = (self.range[d], self.period[d]) {
*x = x.clamp(a.min(b), a.max(b));
}
}
out
}
}
/// Damped Gauss–Newton (Levenberg–Marquardt) for the parameters nearest to `p`, from `start`.
/// Returns the parameters and their distance to `p`.
fn newton(s: &EvalSurface, dom: &Domain, p: Vec3, start: [f64; 2]) -> Option<([f64; 2], f64)> {
let mut uv = dom.clamp(start);
let (mut q, mut su, mut sv) = s.derivs(uv[0], uv[1])?;
let mut dist = q.dist(p);
let mut lambda = 1e-6;
for _ in 0..60 {
let r = p - q;
let (a, b, c) = (su.len2(), su.dot(sv), sv.len2());
let (g0, g1) = (su.dot(r), sv.dot(r));
let scale = a + c;
if !(scale > 0.0) || !scale.is_finite() {
break;
}
let (a2, c2) = (a + lambda * scale, c + lambda * scale);
let det = a2 * c2 - b * b;
if !(det.abs() > 0.0) {
break;
}
let step = [(c2 * g0 - b * g1) / det, (a2 * g1 - b * g0) / det];
let next = dom.clamp([uv[0] + step[0], uv[1] + step[1]]);
let Some((nq, nsu, nsv)) = s.derivs(next[0], next[1]) else { break };
let nd = nq.dist(p);
if nd <= dist {
let moved = (next[0] - uv[0]).abs() + (next[1] - uv[1]).abs();
(uv, q, su, sv, dist) = (next, nq, nsu, nsv, nd);
lambda = (lambda * 0.1).max(1e-12);
if moved <= 1e-15 * (1.0 + uv[0].abs() + uv[1].abs()) {
break;
}
} else {
lambda *= 10.0;
if lambda > 1e8 {
break;
}
}
}
Some((uv, dist))
}
/// Parameters of `p` (a point on the surface within `accept`), from `hint` when that
/// converges, else from the best starts of a grid over the surface's domain.
pub fn project(s: &EvalSurface, dom: &Domain, p: Vec3, hint: Option<[f64; 2]>, accept: f64) -> Option<[f64; 2]> {
if let EvalSurface::Plane { origin, du, dv } = s {
let r = p - *origin;
let (a, b, c) = (du.len2(), du.dot(*dv), dv.len2());
let det = a * c - b * b;
if !(det.abs() > 0.0) {
return None;
}
let (g0, g1) = (du.dot(r), dv.dot(r));
let uv = [(c * g0 - b * g1) / det, (a * g1 - b * g0) / det];
return (uv[0].is_finite() && uv[1].is_finite()).then_some(uv);
}
let mut best: Option<([f64; 2], f64)> = None;
if let Some(h) = hint
&& let Some((uv, d)) = newton(s, dom, p, h)
{
if d <= accept {
return Some(uv);
}
best = Some((uv, d));
}
let (ur, vr) = (dom.range[0]?, dom.range[1]?);
const N: usize = 16;
let mut starts: Vec<(f64, [f64; 2])> = Vec::with_capacity((N + 1) * (N + 1));
for i in 0..=N {
for j in 0..=N {
let uv = [ur.0 + (ur.1 - ur.0) * i as f64 / N as f64, vr.0 + (vr.1 - vr.0) * j as f64 / N as f64];
if let Some((q, _, _)) = s.derivs(uv[0], uv[1]) {
starts.push(((q - p).len2(), uv));
}
}
}
starts.sort_by(|a, b| a.0.total_cmp(&b.0));
for (_, start) in starts.iter().take(4) {
if let Some((uv, d)) = newton(s, dom, p, *start)
&& best.is_none_or(|(_, bd)| d < bd)
{
best = Some((uv, d));
}
}
best.filter(|(_, d)| *d <= accept).map(|(uv, _)| uv)
}
/// `x` moved by whole periods to lie within half a period of `near`.
pub fn unwrap(x: f64, near: f64, period: Option<f64>) -> f64 {
match period {
Some(p) if p > 0.0 => x - ((x - near) / p).round() * p,
_ => x,
}
}
/// Is the surface degenerate along direction `d` at `uv` (a pole: moving along `d` does not
/// move the point)?
pub fn degenerate(s: &EvalSurface, uv: [f64; 2], d: usize) -> bool {
let Some((_, su, sv)) = s.derivs(uv[0], uv[1]) else { return false };
let (this, other) = if d == 0 { (su, sv) } else { (sv, su) };
this.len() <= 1e-9 * other.len().max(1e-12)
}
#[cfg(test)]
mod tests {
use super::*;
use crate::tess::EvalCurve;
fn cylinder() -> EvalSurface {
EvalSurface::Revolved {
profile: EvalCurve { degree: 1, knots: vec![0.0, 0.0, 1.0, 1.0], ctrl: vec![[3.0, 0.0, 0.0, 1.0], [3.0, 0.0, 4.0, 1.0]] },
origin: Vec3::ZERO,
axis: Vec3::Z,
transform: [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]],
swap: false,
}
}
#[test]
fn projects_points_back_to_their_parameters() {
let s = cylinder();
let dom = domain(&s);
for &(u, v) in &[(0.1, 0.2), (0.5, 3.0), (0.9, 6.0), (0.0, 1.0)] {
let (p, _) = s.eval(u, v);
let Some(uv) = project(&s, &dom, p, None, 1e-9) else { panic!("projects") };
let (q, _) = s.eval(uv[0], uv[1]);
assert!(q.dist(p) < 1e-9);
assert!((uv[0] - u).abs() < 1e-9);
}
// Off the surface by more than `accept`: no parameters.
assert!(project(&s, &dom, Vec3::new(5.0, 0.0, 1.0), None, 1e-3).is_none());
}
#[test]
fn unwraps_across_the_period() {
let tau = std::f64::consts::TAU;
assert!((unwrap(0.1, tau - 0.1, Some(tau)) - (tau + 0.1)).abs() < 1e-12);
assert_eq!(unwrap(0.1, 6.0, None), 0.1);
}
}
crates/kernel/src/tess/uv.rs
+2 -187
View File
@@ -2,175 +2,11 @@
//! the plan stores are the ones the evaluators use), and the tolerance-driven divisions of a
//! face's parameter box.
use solidcraft_geom::Vec3;
use solidcraft_geom::tess::EvalSurface;
pub(crate) use solidcraft_geom::tess::uv::{Domain, degenerate, domain, project, unwrap};
use super::size::Rules;
/// The parameter box of a surface and its periods.
#[derive(Clone, Copy, Debug)]
pub(crate) struct Domain {
/// Bounded ranges (`None`: unbounded, as a plane's).
pub range: [Option<(f64, f64)>; 2],
pub period: [Option<f64>; 2],
}
pub(crate) fn domain(s: &EvalSurface) -> Domain {
let tau = std::f64::consts::TAU;
match s {
EvalSurface::Plane { .. } => Domain { range: [None, None], period: [None, None] },
EvalSurface::Nurbs { degree, knots, count, closed, .. } => {
let r = |d: usize| Some((*knots.get(d)?.get(*degree.get(d)?)?, *knots.get(d)?.get(*count.get(d)?)?));
// A direction the surface joins itself along is periodic: a boundary that runs
// across the join is unwrapped past the range's end, and the evaluator wraps it back.
let p = |d: usize| closed.get(d).copied().unwrap_or(false).then(|| r(d).map(|(a, b)| b - a)).flatten().filter(|x| *x > 0.0);
Domain { range: [r(0), r(1)], period: [p(0), p(1)] }
}
EvalSurface::Revolved { profile, swap, .. } => {
let r = (|| Some((*profile.knots.get(profile.degree)?, *profile.knots.get(profile.ctrl.len())?)))();
let (a, b) = ((r, None), (Some((0.0, tau)), Some(tau)));
let (u, v) = if *swap { (b, a) } else { (a, b) };
Domain { range: [u.0, v.0], period: [u.1, v.1] }
}
}
}
impl Domain {
pub(crate) fn clamp(&self, uv: [f64; 2]) -> [f64; 2] {
let mut out = uv;
for (d, x) in out.iter_mut().enumerate() {
if let (Some((a, b)), None) = (self.range[d], self.period[d]) {
*x = x.clamp(a.min(b), a.max(b));
}
}
out
}
}
/// Damped Gauss–Newton (Levenberg–Marquardt) for the parameters nearest to `p`, from `start`.
/// Returns the parameters and their distance to `p`.
fn newton(s: &EvalSurface, dom: &Domain, p: Vec3, start: [f64; 2]) -> Option<([f64; 2], f64)> {
let mut uv = dom.clamp(start);
let (mut q, mut su, mut sv) = s.derivs(uv[0], uv[1])?;
let mut dist = q.dist(p);
let mut lambda = 1e-6;
for _ in 0..60 {
let r = p - q;
let (a, b, c) = (su.len2(), su.dot(sv), sv.len2());
let (g0, g1) = (su.dot(r), sv.dot(r));
let scale = a + c;
if !(scale > 0.0) || !scale.is_finite() {
break;
}
let (a2, c2) = (a + lambda * scale, c + lambda * scale);
let det = a2 * c2 - b * b;
if !(det.abs() > 0.0) {
break;
}
let step = [(c2 * g0 - b * g1) / det, (a2 * g1 - b * g0) / det];
let next = dom.clamp([uv[0] + step[0], uv[1] + step[1]]);
let Some((nq, nsu, nsv)) = s.derivs(next[0], next[1]) else { break };
let nd = nq.dist(p);
if nd <= dist {
let moved = (next[0] - uv[0]).abs() + (next[1] - uv[1]).abs();
(uv, q, su, sv, dist) = (next, nq, nsu, nsv, nd);
lambda = (lambda * 0.1).max(1e-12);
if moved <= 1e-15 * (1.0 + uv[0].abs() + uv[1].abs()) {
break;
}
} else {
lambda *= 10.0;
if lambda > 1e8 {
break;
}
}
}
Some((uv, dist))
}
/// Parameters of `p` (a point on the surface within `accept`), from `hint` when that
/// converges, else from the best starts of a grid over the surface's domain.
pub(crate) fn project(s: &EvalSurface, dom: &Domain, p: Vec3, hint: Option<[f64; 2]>, accept: f64) -> Option<[f64; 2]> {
if let EvalSurface::Plane { origin, du, dv } = s {
let r = p - *origin;
let (a, b, c) = (du.len2(), du.dot(*dv), dv.len2());
let det = a * c - b * b;
if !(det.abs() > 0.0) {
return None;
}
let (g0, g1) = (du.dot(r), dv.dot(r));
let uv = [(c * g0 - b * g1) / det, (a * g1 - b * g0) / det];
return (uv[0].is_finite() && uv[1].is_finite()).then_some(uv);
}
let mut best: Option<([f64; 2], f64)> = None;
if let Some(h) = hint
&& let Some((uv, d)) = newton(s, dom, p, h)
{
if d <= accept {
return Some(uv);
}
best = Some((uv, d));
}
let (ur, vr) = (dom.range[0]?, dom.range[1]?);
let [nu, nv] = start_lines(s);
let mut starts: Vec<(f64, [f64; 2])> = Vec::with_capacity((nu + 1) * (nv + 1));
for i in 0..=nu {
for j in 0..=nv {
let uv = [ur.0 + (ur.1 - ur.0) * i as f64 / nu as f64, vr.0 + (vr.1 - vr.0) * j as f64 / nv as f64];
if let Some((q, _, _)) = s.derivs(uv[0], uv[1]) {
starts.push(((q - p).len2(), uv));
}
}
}
starts.sort_by(|a, b| a.0.total_cmp(&b.0));
for (_, start) in starts.iter().take(4) {
if let Some((uv, d)) = newton(s, dom, p, *start)
&& best.is_none_or(|(_, bd)| d < bd)
{
best = Some((uv, d));
}
}
best.filter(|(_, d)| *d <= accept).map(|(uv, _)| uv)
}
/// Most starts [`project`] tries along one direction, and over the whole domain.
const MAX_START_LINES: usize = 4096;
const MAX_STARTS: usize = 16_384;
/// How many divisions of each direction [`project`] starts from: 16, or two per span of a
/// spline with more. A long spline (a thread's helical flank, as other CAD systems write it to
/// STEP: hundreds of spans along the helix) turns many times between 16 starts, and from the
/// nearest of them Newton settles on another turn, so a point of its boundary was taken as off
/// the surface and the body lost its display plan.
fn start_lines(s: &EvalSurface) -> [usize; 2] {
const N: usize = 16;
let spans = |d: usize| match s {
EvalSurface::Nurbs { degree, count, .. } => count.get(d).zip(degree.get(d)).map_or(0, |(c, g)| c.saturating_sub(*g)),
_ => 0,
};
let mut lines = [0, 1].map(|d| (2 * spans(d)).clamp(N, MAX_START_LINES));
// The longer direction keeps its spans; the other gives way to the cap on starts.
let (long, short) = if lines[0] >= lines[1] { (0, 1) } else { (1, 0) };
lines[short] = lines[short].min((MAX_STARTS / (lines[long] + 1)).max(N));
lines
}
/// `x` moved by whole periods to lie within half a period of `near`.
pub(crate) fn unwrap(x: f64, near: f64, period: Option<f64>) -> f64 {
match period {
Some(p) if p > 0.0 => x - ((x - near) / p).round() * p,
_ => x,
}
}
/// Is the surface degenerate along direction `d` at `uv` (a pole: moving along `d` does not
/// move the point)?
pub(crate) fn degenerate(s: &EvalSurface, uv: [f64; 2], d: usize) -> bool {
let Some((_, su, sv)) = s.derivs(uv[0], uv[1]) else { return false };
let (this, other) = if d == 0 { (su, sv) } else { (sv, su) };
this.len() <= 1e-9 * other.len().max(1e-12)
}
/// Division of `range` along direction `d` so that between neighbouring lines the surface sags
/// at most the tolerance from its chord, the chord is at most `h` long and the normal turns
/// within the rules' facet angle, checked along iso-lines across `across`; at most `max_lines`
@@ -283,6 +119,7 @@ fn limits<'a>(s: &'a EvalSurface, d: usize, across: (f64, f64), rules: &'a Rules
#[cfg(test)]
mod tests {
use super::*;
use solidcraft_geom::Vec3;
use solidcraft_geom::tess::EvalCurve;
fn cylinder() -> EvalSurface {
@@ -295,21 +132,6 @@ mod tests {
}
}
#[test]
fn projects_points_back_to_their_parameters() {
let s = cylinder();
let dom = domain(&s);
for &(u, v) in &[(0.1, 0.2), (0.5, 3.0), (0.9, 6.0), (0.0, 1.0)] {
let (p, _) = s.eval(u, v);
let Some(uv) = project(&s, &dom, p, None, 1e-9) else { panic!("projects") };
let (q, _) = s.eval(uv[0], uv[1]);
assert!(q.dist(p) < 1e-9);
assert!((uv[0] - u).abs() < 1e-9);
}
// Off the surface by more than `accept`: no parameters.
assert!(project(&s, &dom, Vec3::new(5.0, 0.0, 1.0), None, 1e-3).is_none());
}
#[test]
fn divides_the_turn_but_not_the_straight_direction() {
let s = cylinder();
@@ -365,11 +187,4 @@ mod tests {
assert_eq!(lines.first().copied(), Some(a0));
assert_eq!(lines.last().copied(), Some(a1));
}
#[test]
fn unwraps_across_the_period() {
let tau = std::f64::consts::TAU;
assert!((unwrap(0.1, tau - 0.1, Some(tau)) - (tau + 0.1)).abs() < 1e-12);
assert_eq!(unwrap(0.1, 6.0, None), 0.1);
}
}