mirror of
https://github.com/bherbruck/solidcraft.git
synced 2026-11-04 02:55:54 +00:00
brep: pcurves from plain records
solidcraft_geom::brep::pcurve::pcurves(surface, loops, tolerance) computes every edge use's pcurve from our own records: the surface record and, per loop, each use's curve record, interval and direction. No backend's parameter space is involved, so the truck adapter and an SMLib conversion can share it (owner decision 2026-10-21). brep::eval evaluates every schema 4 kind: planes, splines and revolutions through the display evaluator (EvalSurface, EvalCurve, tess::uv), cylinders, cones, spheres and tori in closed form in angle and height, and a surface intersection as the point of both surfaces in the leader's normal plane (what truck's subs converges to), by Newton's method. Exact in our parameters: lines, polylines and whole splines on planes; a circle about the axis of a cylinder, cone, sphere or torus (a u-line); a line along a cylinder's or cone's ruling (a v-line). Each exact pcurve is checked by lifting it; everything else is the cubic B-spline fit the kernel had, with its span and round caps. Curve2::eval moves here from the kernel. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QWhMwdqJYQKQq8RCLSL1LT
This commit is contained in:
jyviko
@@ -10,7 +10,9 @@ use serde::{Deserialize, Serialize};
|
||||
|
||||
mod curve2;
|
||||
pub mod envelope;
|
||||
pub mod eval;
|
||||
pub mod history;
|
||||
pub mod pcurve;
|
||||
mod topology;
|
||||
pub use curve2::Curve2;
|
||||
|
||||
|
||||
@@ -29,6 +29,51 @@ impl Curve2 {
|
||||
let degree = k.len().checked_sub(count.checked_add(1)?)?;
|
||||
Some([*k.get(degree)?, *k.get(count)?])
|
||||
}
|
||||
|
||||
/// The point at `t` (`None` outside the interval or for a malformed record).
|
||||
pub fn eval(&self, t: f64) -> Option<[f64; 2]> {
|
||||
let (knots, h): (&[f64], Vec<[f64; 3]>) = match self {
|
||||
Curve2::Spline { knots, points } => (knots, points.iter().map(|&[u, v]| [u, v, 1.0]).collect()),
|
||||
Curve2::RationalSpline { knots, homogeneous } => (knots, homogeneous.clone()),
|
||||
};
|
||||
let p = knots.len().checked_sub(h.len() + 1)?;
|
||||
let span = crate::bspline::find_span(knots, p, h.len(), t)?;
|
||||
let n = basis(knots, p, span, t)?;
|
||||
let mut s = [0.0; 3];
|
||||
for (j, nj) in n.iter().enumerate() {
|
||||
let q = h.get(span + j - p)?;
|
||||
for (a, x) in s.iter_mut().zip(q) {
|
||||
*a += nj * x;
|
||||
}
|
||||
}
|
||||
(s[2].abs() > 0.0).then(|| [s[0] / s[2], s[1] / s[2]])
|
||||
}
|
||||
}
|
||||
|
||||
/// The `p + 1` basis functions non-zero on knot span `span` at `t`, of any degree (NURBS Book
|
||||
/// A2.2).
|
||||
fn basis(knots: &[f64], p: usize, span: usize, t: f64) -> Option<Vec<f64>> {
|
||||
if span < p || span + p >= knots.len() {
|
||||
return None;
|
||||
}
|
||||
let k = |i: usize| knots.get(i).copied().unwrap_or(0.0);
|
||||
let mut n = vec![0.0; p + 1];
|
||||
let (mut left, mut right) = (vec![0.0; p + 1], vec![0.0; p + 1]);
|
||||
*n.first_mut()? = 1.0;
|
||||
for j in 1..=p {
|
||||
*left.get_mut(j)? = t - k(span + 1 - j);
|
||||
*right.get_mut(j)? = k(span + j) - t;
|
||||
let mut saved = 0.0;
|
||||
for r in 0..j {
|
||||
let (rr, lj) = (*right.get(r + 1)?, *left.get(j - r)?);
|
||||
let den = rr + lj;
|
||||
let tmp = if den != 0.0 { n.get(r)? / den } else { 0.0 };
|
||||
*n.get_mut(r)? = saved + rr * tmp;
|
||||
saved = lj * tmp;
|
||||
}
|
||||
*n.get_mut(j)? = saved;
|
||||
}
|
||||
Some(n)
|
||||
}
|
||||
|
||||
pub(super) fn curve2(c: &Curve2, budget: &mut Budget) -> Result<()> {
|
||||
|
||||
@@ -0,0 +1,390 @@
|
||||
//! Points of curve and surface records, in the records' own parameters (docs/native-brep.md,
|
||||
//! "Geometry"), and the parameters of a point on a surface. Planes, splines and revolutions go
|
||||
//! through the display evaluator ([`EvalSurface`], [`EvalCurve`], [`uv`]), which parameterises
|
||||
//! them as their producer does; cylinders, cones, spheres and tori are evaluated and inverted in
|
||||
//! closed form, in angle and height. A `surface_intersection` is the point of both surfaces in
|
||||
//! the plane through its leader's point normal to the leader there, which truck's `subs`
|
||||
//! converges to (within its 10⁻⁶ mm), found by Newton's method on the three planes.
|
||||
//!
|
||||
//! Call on validated records. Evaluation of a well-formed record never panics: a point that
|
||||
//! cannot be had (a vanishing weight, a Newton step that does not converge) is `None`.
|
||||
use super::pcurve::Error;
|
||||
use super::{Curve, Frame, Surface};
|
||||
use crate::Vec3;
|
||||
use crate::bspline::MAX_SPLINE_DEGREE;
|
||||
use crate::tess::uv::{self, Domain};
|
||||
use crate::tess::{EvalCurve, EvalSurface, MAX_EVAL_KNOTS, MAX_EVAL_POINTS};
|
||||
use std::cell::Cell;
|
||||
use std::f64::consts::{FRAC_PI_2, TAU};
|
||||
|
||||
/// Newton steps of one intersection point.
|
||||
const MAX_NEWTON: usize = 64;
|
||||
|
||||
fn v3(p: [f64; 3]) -> Vec3 {
|
||||
Vec3::new(p[0], p[1], p[2])
|
||||
}
|
||||
fn unsupported(what: impl std::fmt::Display) -> Error {
|
||||
Error::Failed(format!("not supported yet: {what}"))
|
||||
}
|
||||
|
||||
/// A frame's origin and axes, z = x × y.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub struct Axes {
|
||||
pub(crate) o: Vec3,
|
||||
pub(crate) x: Vec3,
|
||||
pub(crate) y: Vec3,
|
||||
pub(crate) z: Vec3,
|
||||
}
|
||||
impl Axes {
|
||||
pub(crate) fn new(f: &Frame) -> Self {
|
||||
let (x, y) = (v3(f.x_axis), v3(f.y_axis));
|
||||
Axes { o: v3(f.origin), x, y, z: x.cross(y) }
|
||||
}
|
||||
/// `p` in the frame: (x, y, z) coordinates.
|
||||
pub(crate) fn local(&self, p: Vec3) -> [f64; 3] {
|
||||
let q = p - self.o;
|
||||
[q.dot(self.x), q.dot(self.y), q.dot(self.z)]
|
||||
}
|
||||
/// The radial direction at angle `u` and its derivative.
|
||||
fn radial(&self, u: f64) -> (Vec3, Vec3) {
|
||||
let (s, c) = u.sin_cos();
|
||||
(self.x * c + self.y * s, self.y * c - self.x * s)
|
||||
}
|
||||
}
|
||||
|
||||
/// The angle of (a, b) in [0, 2π), `near`'s (or 0) on the axis, where it is undefined.
|
||||
fn angle(a: f64, b: f64, near: Option<f64>) -> f64 {
|
||||
if a == 0.0 && b == 0.0 {
|
||||
return near.unwrap_or(0.0);
|
||||
}
|
||||
let u = b.atan2(a);
|
||||
if u < 0.0 { u + TAU } else { u }
|
||||
}
|
||||
|
||||
/// A surface record ready to evaluate.
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
pub enum SurfaceEval {
|
||||
/// A plane, spline or revolution, on the display evaluator.
|
||||
Eval(EvalSurface, Domain),
|
||||
Cylinder {
|
||||
axes: Axes,
|
||||
radius: f64,
|
||||
},
|
||||
Cone {
|
||||
axes: Axes,
|
||||
radius_at_origin: f64,
|
||||
slope: f64,
|
||||
},
|
||||
Sphere {
|
||||
axes: Axes,
|
||||
radius: f64,
|
||||
},
|
||||
Torus {
|
||||
axes: Axes,
|
||||
major: f64,
|
||||
minor: f64,
|
||||
},
|
||||
}
|
||||
|
||||
impl SurfaceEval {
|
||||
pub fn new(s: &Surface) -> Result<Self, Error> {
|
||||
let eval = |e: EvalSurface| {
|
||||
let dom = uv::domain(&e);
|
||||
Ok(SurfaceEval::Eval(e, dom))
|
||||
};
|
||||
match s {
|
||||
Surface::Plane { origin, u_axis, v_axis } => eval(EvalSurface::Plane { origin: v3(*origin), du: v3(*u_axis), dv: v3(*v_axis) }),
|
||||
Surface::PlaneThrough { origin, u_point, v_point } => {
|
||||
let o = v3(*origin);
|
||||
eval(EvalSurface::Plane { origin: o, du: v3(*u_point) - o, dv: v3(*v_point) - o })
|
||||
}
|
||||
Surface::Cylinder { frame, radius } => Ok(SurfaceEval::Cylinder { axes: Axes::new(frame), radius: *radius }),
|
||||
Surface::Cone { frame, radius_at_origin, half_angle } => {
|
||||
Ok(SurfaceEval::Cone { axes: Axes::new(frame), radius_at_origin: *radius_at_origin, slope: half_angle.tan() })
|
||||
}
|
||||
Surface::Sphere { frame, radius } => Ok(SurfaceEval::Sphere { axes: Axes::new(frame), radius: *radius }),
|
||||
Surface::Torus { frame, major_radius, minor_radius } => {
|
||||
Ok(SurfaceEval::Torus { axes: Axes::new(frame), major: *major_radius, minor: *minor_radius })
|
||||
}
|
||||
Surface::Spline { u_knots, v_knots, points } => {
|
||||
let rows: Vec<Vec<[f64; 4]>> = points.iter().map(|r| r.iter().map(|p| [p[0], p[1], p[2], 1.0]).collect()).collect();
|
||||
eval(nurbs(u_knots, v_knots, &rows)?)
|
||||
}
|
||||
Surface::RationalSpline { u_knots, v_knots, homogeneous } => eval(nurbs(u_knots, v_knots, homogeneous)?),
|
||||
Surface::PlacedRevolution { profile, origin, axis, matrix, forward } => {
|
||||
let profile = match CurveEval::new(profile)? {
|
||||
CurveEval::Spline(c) => c,
|
||||
CurveEval::Line { start, end } => EvalCurve {
|
||||
degree: 1,
|
||||
knots: vec![0.0, 0.0, 1.0, 1.0],
|
||||
ctrl: vec![[start.x, start.y, start.z, 1.0], [end.x, end.y, end.z, 1.0]],
|
||||
},
|
||||
_ => return Err(unsupported("a revolution of a conic, polyline or intersection profile")),
|
||||
};
|
||||
let axis = v3(*axis).normalized().ok_or_else(|| Error::Invalid("revolution axis".into()))?;
|
||||
// Column-major: matrix[column][row].
|
||||
let transform = [0, 1, 2].map(|r| [matrix[0][r], matrix[1][r], matrix[2][r], matrix[3][r]]);
|
||||
eval(EvalSurface::Revolved { profile, origin: v3(*origin), axis, transform, swap: !forward })
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// The parameter box and periods.
|
||||
pub fn domain(&self) -> Domain {
|
||||
let turn = (Some((0.0, TAU)), Some(TAU));
|
||||
let open = (None, None);
|
||||
let ((ru, pu), (rv, pv)) = match self {
|
||||
SurfaceEval::Eval(_, dom) => return *dom,
|
||||
SurfaceEval::Cylinder { .. } | SurfaceEval::Cone { .. } => (turn, open),
|
||||
SurfaceEval::Sphere { .. } => (turn, (Some((-FRAC_PI_2, FRAC_PI_2)), None)),
|
||||
SurfaceEval::Torus { .. } => (turn, turn),
|
||||
};
|
||||
Domain { range: [ru, rv], period: [pu, pv] }
|
||||
}
|
||||
|
||||
/// Point and first partial derivatives at (u, v).
|
||||
pub fn derivs(&self, u: f64, v: f64) -> Option<(Vec3, Vec3, Vec3)> {
|
||||
if !(u.is_finite() && v.is_finite()) {
|
||||
return None;
|
||||
}
|
||||
Some(match self {
|
||||
SurfaceEval::Eval(e, _) => return e.derivs(u, v),
|
||||
SurfaceEval::Cylinder { axes: a, radius } => {
|
||||
let (r, dr) = a.radial(u);
|
||||
(a.o + a.z * v + r * *radius, dr * *radius, a.z)
|
||||
}
|
||||
SurfaceEval::Cone { axes: a, radius_at_origin, slope } => {
|
||||
let (r, dr) = a.radial(u);
|
||||
let rv = radius_at_origin + v * slope;
|
||||
(a.o + a.z * v + r * rv, dr * rv, a.z + r * *slope)
|
||||
}
|
||||
SurfaceEval::Sphere { axes: a, radius } => {
|
||||
let (r, dr) = a.radial(u);
|
||||
let (s, c) = v.sin_cos();
|
||||
(a.o + (r * c + a.z * s) * *radius, dr * (c * radius), (a.z * c - r * s) * *radius)
|
||||
}
|
||||
SurfaceEval::Torus { axes: a, major, minor } => {
|
||||
let (r, dr) = a.radial(u);
|
||||
let (s, c) = v.sin_cos();
|
||||
let ring = major + minor * c;
|
||||
(a.o + r * ring + a.z * (minor * s), dr * ring, (a.z * c - r * s) * *minor)
|
||||
}
|
||||
})
|
||||
}
|
||||
|
||||
pub fn point(&self, uv: [f64; 2]) -> Option<Vec3> {
|
||||
self.derivs(uv[0], uv[1]).map(|d| d.0)
|
||||
}
|
||||
|
||||
/// The parameters of `p`, a point within `accept` of the surface, from `near` when given
|
||||
/// (not unwrapped to it).
|
||||
pub fn project(&self, p: Vec3, near: Option<[f64; 2]>, accept: f64) -> Option<[f64; 2]> {
|
||||
let q = match self {
|
||||
SurfaceEval::Eval(e, dom) => return uv::project(e, dom, p, near, accept),
|
||||
SurfaceEval::Cylinder { axes, .. } => {
|
||||
let [a, b, h] = axes.local(p);
|
||||
[angle(a, b, near.map(|n| n[0])), h]
|
||||
}
|
||||
SurfaceEval::Cone { axes, radius_at_origin, slope } => {
|
||||
// The nearest point of the generator line through p's half-plane.
|
||||
let [a, b, h] = axes.local(p);
|
||||
[angle(a, b, near.map(|n| n[0])), ((a.hypot(b) - radius_at_origin) * slope + h) / (1.0 + slope * slope)]
|
||||
}
|
||||
SurfaceEval::Sphere { axes, .. } => {
|
||||
let [a, b, h] = axes.local(p);
|
||||
[angle(a, b, near.map(|n| n[0])), h.atan2(a.hypot(b))]
|
||||
}
|
||||
SurfaceEval::Torus { axes, major, .. } => {
|
||||
let [a, b, h] = axes.local(p);
|
||||
let v = h.atan2(a.hypot(b) - major);
|
||||
[angle(a, b, near.map(|n| n[0])), if v < 0.0 { v + TAU } else { v }]
|
||||
}
|
||||
};
|
||||
(q[0].is_finite() && q[1].is_finite() && self.point(q).is_some_and(|s| s.dist(p) <= accept)).then_some(q)
|
||||
}
|
||||
|
||||
/// Whether moving along each direction at `uv` leaves the point where it is (a pole).
|
||||
pub fn pole(&self, uv: [f64; 2]) -> [bool; 2] {
|
||||
let Some((_, su, sv)) = self.derivs(uv[0], uv[1]) else { return [false; 2] };
|
||||
let (a, b) = (su.len(), sv.len());
|
||||
[a <= 1e-9 * b.max(1e-12), b <= 1e-9 * a.max(1e-12)]
|
||||
}
|
||||
}
|
||||
|
||||
/// A tensor spline's rows (`rows[i][j]`, i along u) as the evaluator's u-fastest array.
|
||||
fn nurbs(u_knots: &[f64], v_knots: &[f64], rows: &[Vec<[f64; 4]>]) -> Result<EvalSurface, Error> {
|
||||
let (nu, nv) = (rows.len(), rows.first().map_or(0, Vec::len));
|
||||
let degree = |k: &[f64], n: usize| k.len().checked_sub(n + 1).ok_or_else(|| Error::Invalid("spline degree".into()));
|
||||
let (pu, pv) = (degree(u_knots, nu)?, degree(v_knots, nv)?);
|
||||
if pu.max(pv) > MAX_SPLINE_DEGREE || pu.min(pv) == 0 {
|
||||
return Err(unsupported(format!("a spline surface of degree ({pu}, {pv})")));
|
||||
}
|
||||
if u_knots.len().max(v_knots.len()) > MAX_EVAL_KNOTS || nu.saturating_mul(nv) > MAX_EVAL_POINTS {
|
||||
return Err(Error::TooLarge(format!("a spline surface of {nu} × {nv} points")));
|
||||
}
|
||||
let mut ctrl = Vec::with_capacity(nu * nv);
|
||||
for j in 0..nv {
|
||||
for row in rows {
|
||||
ctrl.push(*row.get(j).ok_or_else(|| Error::Invalid("ragged surface".into()))?);
|
||||
}
|
||||
}
|
||||
Ok(EvalSurface::Nurbs { degree: [pu, pv], knots: [u_knots.to_vec(), v_knots.to_vec()], count: [nu, nv], ctrl })
|
||||
}
|
||||
|
||||
/// A curve record ready to evaluate.
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
pub enum CurveEval {
|
||||
Line {
|
||||
start: Vec3,
|
||||
end: Vec3,
|
||||
},
|
||||
/// A circle or an ellipse.
|
||||
Conic {
|
||||
axes: Axes,
|
||||
radii: [f64; 2],
|
||||
},
|
||||
Polyline(Vec<Vec3>),
|
||||
Spline(EvalCurve),
|
||||
Intersection(Box<Intersection>),
|
||||
}
|
||||
|
||||
/// A `surface_intersection` ready to evaluate.
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
pub struct Intersection {
|
||||
surfaces: [SurfaceEval; 2],
|
||||
leader: CurveEval,
|
||||
tolerance: f64,
|
||||
/// Both surfaces' parameters of the last point found: the next point starts from them.
|
||||
hint: Cell<Option<[[f64; 2]; 2]>>,
|
||||
}
|
||||
|
||||
impl CurveEval {
|
||||
pub fn new(c: &Curve) -> Result<Self, Error> {
|
||||
Ok(match c {
|
||||
Curve::Line { start, end } => CurveEval::Line { start: v3(*start), end: v3(*end) },
|
||||
Curve::Circle { frame, radius } => CurveEval::Conic { axes: Axes::new(frame), radii: [*radius; 2] },
|
||||
Curve::Ellipse { frame, radii } => CurveEval::Conic { axes: Axes::new(frame), radii: *radii },
|
||||
Curve::Polyline { points } => CurveEval::Polyline(points.iter().copied().map(v3).collect()),
|
||||
Curve::Spline { knots, points } => spline(knots, points.iter().map(|p| [p[0], p[1], p[2], 1.0]).collect())?,
|
||||
Curve::RationalSpline { knots, homogeneous } => spline(knots, homogeneous.clone())?,
|
||||
Curve::SurfaceIntersection { surfaces, leader, tolerance } => {
|
||||
let leader = CurveEval::new(leader)?;
|
||||
if !matches!(leader, CurveEval::Polyline(_) | CurveEval::Spline(_)) {
|
||||
return Err(Error::Invalid("intersection leader".into()));
|
||||
}
|
||||
CurveEval::Intersection(Box::new(Intersection {
|
||||
surfaces: [SurfaceEval::new(&surfaces[0])?, SurfaceEval::new(&surfaces[1])?],
|
||||
leader,
|
||||
tolerance: *tolerance,
|
||||
hint: Cell::new(None),
|
||||
}))
|
||||
}
|
||||
})
|
||||
}
|
||||
|
||||
/// The point at `t`.
|
||||
pub fn point(&self, t: f64) -> Option<Vec3> {
|
||||
if !t.is_finite() {
|
||||
return None;
|
||||
}
|
||||
match self {
|
||||
CurveEval::Line { start, end } => Some(*start + (*end - *start) * t),
|
||||
CurveEval::Conic { axes: a, radii } => {
|
||||
let (s, c) = t.sin_cos();
|
||||
Some(a.o + a.x * (radii[0] * c) + a.y * (radii[1] * s))
|
||||
}
|
||||
CurveEval::Polyline(_) | CurveEval::Spline(_) => self.leader(t).map(|d| d.0),
|
||||
CurveEval::Intersection(x) => x.point(t),
|
||||
}
|
||||
}
|
||||
|
||||
/// Point and derivative of a polyline (truck's: the segment from ⌊t⌋, the last one at the
|
||||
/// end) or a spline.
|
||||
fn leader(&self, t: f64) -> Option<(Vec3, Vec3)> {
|
||||
match self {
|
||||
CurveEval::Polyline(p) => {
|
||||
let n = p.len();
|
||||
let last = *p.last()?;
|
||||
if t <= 0.0 {
|
||||
return Some((*p.first()?, *p.get(1)? - *p.first()?));
|
||||
}
|
||||
if t + 1.0 >= n as f64 {
|
||||
return Some((last, last - *p.get(n.checked_sub(2)?)?));
|
||||
}
|
||||
let i = t as usize;
|
||||
let (a, b) = (*p.get(i)?, *p.get(i + 1)?);
|
||||
Some((a + (b - a) * (t - i as f64), b - a))
|
||||
}
|
||||
CurveEval::Spline(c) => c.eval(t),
|
||||
_ => None,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn spline(knots: &[f64], ctrl: Vec<[f64; 4]>) -> Result<CurveEval, Error> {
|
||||
let degree = knots.len().checked_sub(ctrl.len() + 1).ok_or_else(|| Error::Invalid("spline degree".into()))?;
|
||||
if !(1..=MAX_SPLINE_DEGREE).contains(°ree) {
|
||||
return Err(unsupported(format!("a spline curve of degree {degree}")));
|
||||
}
|
||||
Ok(CurveEval::Spline(EvalCurve { degree, knots: knots.to_vec(), ctrl }))
|
||||
}
|
||||
|
||||
impl Intersection {
|
||||
/// The point in the leader's normal plane at `t`: from the last point's parameters when that
|
||||
/// lands within the trace tolerance of the leader (as truck takes a hinted point), else from
|
||||
/// the leader's point alone.
|
||||
fn point(&self, t: f64) -> Option<Vec3> {
|
||||
let (lead, along) = self.leader.leader(t)?;
|
||||
if !(along.len2() > 0.0) {
|
||||
return None;
|
||||
}
|
||||
let near = self.tolerance.max(1e-6);
|
||||
let hinted = self.hint.get().and_then(|h| self.solve(lead, along, Some(h))).filter(|(x, _)| x.dist(lead) <= near);
|
||||
let (x, uv) = match hinted {
|
||||
Some(found) => found,
|
||||
None => self.solve(lead, along, None)?,
|
||||
};
|
||||
self.hint.set(Some(uv));
|
||||
Some(x)
|
||||
}
|
||||
|
||||
/// Newton on the planes: each surface's tangent plane at the foot of the current point, and
|
||||
/// the leader's normal plane, until the point stops moving.
|
||||
fn solve(&self, lead: Vec3, along: Vec3, hint: Option<[[f64; 2]; 2]>) -> Option<(Vec3, [[f64; 2]; 2])> {
|
||||
let mut x = lead;
|
||||
let mut uv = hint;
|
||||
for _ in 0..MAX_NEWTON {
|
||||
let mut planes = [(Vec3::ZERO, 0.0); 2];
|
||||
let mut feet = [[0.0; 2]; 2];
|
||||
for (i, s) in self.surfaces.iter().enumerate() {
|
||||
let q = s.project(x, uv.and_then(|h| h.get(i).copied()), f64::INFINITY)?;
|
||||
let (f, su, sv) = s.derivs(q[0], q[1])?;
|
||||
let n = su.cross(sv).normalized()?;
|
||||
*planes.get_mut(i)? = (n, n.dot(f));
|
||||
*feet.get_mut(i)? = q;
|
||||
}
|
||||
uv = Some(feet);
|
||||
let [(n0, d0), (n1, d1)] = planes;
|
||||
let y = solve3([n0, n1, along], [d0, d1, along.dot(lead)])?;
|
||||
let step = y.dist(x);
|
||||
x = y;
|
||||
if step <= 1e-12 * (1.0 + x.len()) {
|
||||
return Some((x, feet));
|
||||
}
|
||||
}
|
||||
None
|
||||
}
|
||||
}
|
||||
|
||||
/// The point on the three planes `rows[i] · x = rhs[i]` (Cramer's rule), `None` when they meet
|
||||
/// in no single point (tangent surfaces).
|
||||
fn solve3(rows: [Vec3; 3], rhs: [f64; 3]) -> Option<Vec3> {
|
||||
let [a, b, c] = rows;
|
||||
let (bc, ca, ab) = (b.cross(c), c.cross(a), a.cross(b));
|
||||
let det = a.dot(bc);
|
||||
if !(det.abs() > 1e-12 * a.len() * b.len() * c.len()) {
|
||||
return None;
|
||||
}
|
||||
let x = (bc * rhs[0] + ca * rhs[1] + ab * rhs[2]) * (1.0 / det);
|
||||
x.is_finite().then_some(x)
|
||||
}
|
||||
@@ -0,0 +1,621 @@
|
||||
//! Pcurves from plain records: each edge use's curve in its face's surface parameters, on the
|
||||
//! edge's own interval and running in its curve's direction (docs/native-brep.md, "Pcurves"),
|
||||
//! within a tolerance when lifted onto the surface. One routine for every backend: computed
|
||||
//! from the 3D edge and our surface record, never taken from a kernel's own parameter space.
|
||||
//!
|
||||
//! Exact where the edge lies on its surface in closed form: a line, polyline or spline on a
|
||||
//! plane (the plane's parameters are affine in space: the curve's knots and weights with its
|
||||
//! control points mapped); a circle about the axis of a cylinder, cone, sphere or torus (v
|
||||
//! constant, u the circle's angle plus a constant, or minus it where the circle turns the other
|
||||
//! way); a line along a ruling of a cylinder or cone (u constant, v linear). Each exact pcurve is
|
||||
//! checked by lifting it; anything else, and any exact one that misses, is a cubic B-spline fitted
|
||||
//! within the tolerance, never a polyline.
|
||||
//!
|
||||
//! Periodic parameters are unwrapped along each curve, and each loop's uses are moved by whole
|
||||
//! periods so every use starts where the one before it ends: a seam edge used twice gets two
|
||||
//! pcurves a period apart. At a pole the parameter along the collapsed direction is undefined;
|
||||
//! there it is taken from the nearest sample where it is defined, so a pcurve runs straight into
|
||||
//! the pole, and where two uses meet at a pole the loop runs along the pole's line.
|
||||
use super::eval::{CurveEval, SurfaceEval};
|
||||
use super::{Budget, Curve, Curve2, Surface, curve, curve_natural, interval, surface};
|
||||
use crate::Vec3;
|
||||
use crate::tess::uv::{self, Domain};
|
||||
|
||||
/// Why no pcurves: the records are malformed or out of bounds (`Invalid`), no pcurve within
|
||||
/// the tolerance was found or the geometry is not supported (`Failed`), or the fit would exceed
|
||||
/// its caps (`TooLarge`).
|
||||
#[derive(Clone, Debug, PartialEq, thiserror::Error)]
|
||||
pub enum Error {
|
||||
#[error("{0}")]
|
||||
Invalid(String),
|
||||
#[error("{0}")]
|
||||
Failed(String),
|
||||
#[error("{0}")]
|
||||
TooLarge(String),
|
||||
}
|
||||
|
||||
type Result<T> = std::result::Result<T, Error>;
|
||||
|
||||
/// One oriented edge of a loop, as [`pcurves`] reads it: the edge's index (a seam edge, used
|
||||
/// twice, is fitted once), its curve record and interval, and whether the loop runs along it.
|
||||
#[derive(Clone, Copy, Debug)]
|
||||
pub struct EdgeUse<'a> {
|
||||
pub edge: u32,
|
||||
pub curve: &'a Curve,
|
||||
pub interval: [f64; 2],
|
||||
pub forward: bool,
|
||||
}
|
||||
|
||||
/// 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.
|
||||
const CHECKS: usize = 16;
|
||||
|
||||
/// For each loop of a face on `surface`, the pcurve of each of its edge uses, in order, within
|
||||
/// `tolerance` mm of the edge when lifted onto the surface. Loops are as the record keeps them
|
||||
/// (relative to the surface's natural normal: outer counter-clockwise, inner clockwise); each
|
||||
/// pcurve runs on its use's interval in its curve's direction, whatever `forward` says.
|
||||
pub fn pcurves(surface_record: &Surface, loops: &[Vec<EdgeUse<'_>>], tolerance: f64) -> Result<Vec<Vec<Curve2>>> {
|
||||
if !(tolerance.is_finite() && tolerance > 0.0 && tolerance <= 1e3) {
|
||||
return Err(Error::Invalid(format!("pcurve tolerance {tolerance}")));
|
||||
}
|
||||
let uses = loops.iter().map(Vec::len).fold(0usize, usize::saturating_add);
|
||||
if uses > MAX_USES {
|
||||
return Err(Error::TooLarge(format!("a face with {uses} edge uses")));
|
||||
}
|
||||
let mut budget = Budget(0);
|
||||
surface(surface_record, 0, &mut budget).map_err(|e| Error::Invalid(e.to_string()))?;
|
||||
let on = On::new(surface_record)?;
|
||||
// A seam edge used twice is fitted once.
|
||||
let mut done: Vec<(u32, Curve2)> = Vec::new();
|
||||
let mut out = Vec::with_capacity(loops.len());
|
||||
for (w, l) in loops.iter().enumerate() {
|
||||
let mut fitted = Vec::with_capacity(l.len());
|
||||
for (i, u) in l.iter().enumerate() {
|
||||
let pc = match done.iter().find(|(e, _)| *e == u.edge) {
|
||||
Some((_, pc)) => pc.clone(),
|
||||
None => {
|
||||
curve(u.curve, 0, &mut budget).map_err(|e| Error::Invalid(format!("loop {w} use {i}: {e}")))?;
|
||||
interval(&u.interval, curve_natural(u.curve)).map_err(|e| Error::Invalid(format!("loop {w} use {i}: {e}")))?;
|
||||
let c = CurveEval::new(u.curve)?;
|
||||
let pc = on.pcurve(u.curve, &c, u.interval, tolerance).map_err(|e| {
|
||||
let at = |t: f64| c.point(t).map(|p| [p.x, p.y, p.z]);
|
||||
let message = format!("pcurve of loop {w} use {i} ({:?} to {:?}): {e}", at(u.interval[0]), at(u.interval[1]));
|
||||
match e {
|
||||
Error::Invalid(_) => Error::Invalid(message),
|
||||
Error::Failed(_) => Error::Failed(message),
|
||||
Error::TooLarge(_) => Error::TooLarge(message),
|
||||
}
|
||||
})?;
|
||||
done.push((u.edge, pc.clone()));
|
||||
pc
|
||||
}
|
||||
};
|
||||
fitted.push((pc, u.forward));
|
||||
}
|
||||
on.close(&mut fitted)?;
|
||||
out.push(fitted.into_iter().map(|(c, _)| c).collect());
|
||||
}
|
||||
Ok(out)
|
||||
}
|
||||
|
||||
/// A face's surface: its evaluator and parameter box.
|
||||
struct On {
|
||||
surface: SurfaceEval,
|
||||
dom: Domain,
|
||||
}
|
||||
|
||||
/// 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)?;
|
||||
let dom = surface.domain();
|
||||
Ok(On { surface, dom })
|
||||
}
|
||||
|
||||
fn pcurve(&self, record: &Curve, curve: &CurveEval, range: [f64; 2], tol: f64) -> Result<Curve2> {
|
||||
if let Some(pc) = self.exact(record, curve, range)?
|
||||
&& self.lifts(&pc, curve, range, 0.5 * tol)
|
||||
{
|
||||
return Ok(pc);
|
||||
}
|
||||
self.fit(record, curve, range, tol)
|
||||
}
|
||||
|
||||
/// The closed-form pcurve where there is one (unchecked).
|
||||
fn exact(&self, record: &Curve, curve: &CurveEval, [t0, t1]: [f64; 2]) -> Result<Option<Curve2>> {
|
||||
let segment = |a: [f64; 2], b: [f64; 2]| Curve2::Spline { knots: vec![t0, t0, t1, t1], points: vec![a, b] };
|
||||
let at = |t: f64| curve.point(t).ok_or_else(|| Error::Failed(format!("the curve has no point at {t}")));
|
||||
Ok(match (&self.surface, record) {
|
||||
(SurfaceEval::Eval(plane @ crate::tess::EvalSurface::Plane { .. }, dom), _) => {
|
||||
let map = |p: Vec3| uv::project(plane, dom, p, None, 0.0).ok_or_else(|| Error::Invalid("degenerate plane".into()));
|
||||
match record {
|
||||
// The line's own ends where the interval reaches them.
|
||||
Curve::Line { start, end } => {
|
||||
let end_at = |t: f64| {
|
||||
if t == 0.0 {
|
||||
Ok(Vec3::from(*start))
|
||||
} else if t == 1.0 {
|
||||
Ok(Vec3::from(*end))
|
||||
} else {
|
||||
at(t)
|
||||
}
|
||||
};
|
||||
Some(segment(map(end_at(t0)?)?, map(end_at(t1)?)?))
|
||||
}
|
||||
Curve::Polyline { .. } => {
|
||||
let mut ts = vec![t0];
|
||||
ts.extend((t0.floor() as i64 + 1..).map(|k| k as f64).take_while(|k| *k < t1).take(MAX_SPANS));
|
||||
ts.push(t1);
|
||||
if ts.len() > 4096 {
|
||||
return Ok(None);
|
||||
}
|
||||
let mut knots = vec![t0];
|
||||
knots.extend(ts.iter().copied());
|
||||
knots.push(t1);
|
||||
Some(Curve2::Spline { knots, points: ts.iter().map(|t| map(at(*t)?)).collect::<Result<_>>()? })
|
||||
}
|
||||
Curve::Spline { knots, points } if curve_natural(record) == super::Natural::Range(t0, t1) => {
|
||||
Some(Curve2::Spline { knots: knots.clone(), points: points.iter().map(|p| map(Vec3::from(*p))).collect::<Result<_>>()? })
|
||||
}
|
||||
Curve::RationalSpline { knots, homogeneous } if curve_natural(record) == super::Natural::Range(t0, t1) => {
|
||||
Some(Curve2::RationalSpline {
|
||||
knots: knots.clone(),
|
||||
homogeneous: homogeneous
|
||||
.iter()
|
||||
.map(|h| {
|
||||
let [u, v] = map(Vec3::new(h[0] / h[3], h[1] / h[3], h[2] / h[3]))?;
|
||||
Ok([u * h[3], v * h[3], h[3]])
|
||||
})
|
||||
.collect::<Result<_>>()?,
|
||||
})
|
||||
}
|
||||
_ => None,
|
||||
}
|
||||
}
|
||||
(
|
||||
SurfaceEval::Cylinder { axes, .. }
|
||||
| SurfaceEval::Cone { axes, .. }
|
||||
| SurfaceEval::Sphere { axes, .. }
|
||||
| SurfaceEval::Torus { axes, .. },
|
||||
Curve::Circle { frame, radius },
|
||||
) => {
|
||||
// About the axis: u the circle's angle from where its x axis points, v fixed.
|
||||
let c = super::eval::Axes::new(frame);
|
||||
let turn = c.z.dot(axes.z);
|
||||
if (turn.abs() - 1.0).abs() > 1e-9 {
|
||||
return Ok(None);
|
||||
}
|
||||
let sign = turn.signum();
|
||||
let u0 = c.x.dot(axes.y).atan2(c.x.dot(axes.x));
|
||||
let h = axes.local(c.o)[2];
|
||||
let v = match &self.surface {
|
||||
SurfaceEval::Sphere { .. } => h.atan2(*radius),
|
||||
SurfaceEval::Torus { major, .. } => {
|
||||
let v = h.atan2(radius - major);
|
||||
if v < 0.0 { v + std::f64::consts::TAU } else { v }
|
||||
}
|
||||
_ => h,
|
||||
};
|
||||
Some(segment([u0 + sign * t0, v], [u0 + sign * t1, v]))
|
||||
}
|
||||
(SurfaceEval::Cylinder { axes, .. } | SurfaceEval::Cone { axes, .. }, Curve::Line { .. }) => {
|
||||
// Along a ruling: u where the line is, v its height.
|
||||
let [a, b, _] = axes.local(at(0.5 * (t0 + t1))?);
|
||||
let height = |t: f64| at(t).map(|p| axes.local(p)[2]);
|
||||
let u = b.atan2(a);
|
||||
Some(segment([u, height(t0)?], [u, height(t1)?]))
|
||||
}
|
||||
_ => None,
|
||||
})
|
||||
}
|
||||
|
||||
/// Whether `pc` lifts onto `curve` within `target` along `range` (evenly, and inside its
|
||||
/// spans).
|
||||
fn lifts(&self, pc: &Curve2, curve: &CurveEval, [t0, t1]: [f64; 2], target: f64) -> bool {
|
||||
let knots = match pc {
|
||||
Curve2::Spline { knots, .. } | Curve2::RationalSpline { knots, .. } => knots,
|
||||
};
|
||||
let even = (0..=CHECKS).map(|k| t0 + (t1 - t0) * k as f64 / CHECKS as f64);
|
||||
let spans = knots.windows(2).filter_map(|k| Some((*k.first()?, *k.get(1)?))).filter(|(a, b)| a < b).take(1 << 12).map(|(a, b)| 0.5 * (a + b));
|
||||
even.chain(spans).all(|t| {
|
||||
let lifted = pc.eval(t).and_then(|q| self.surface.point(q));
|
||||
matches!((lifted, curve.point(t)), (Some(a), Some(b)) if a.dist(b) <= target)
|
||||
})
|
||||
}
|
||||
|
||||
/// The parameters of `p` within `accept`, from `near` (and unwrapped to it) when given.
|
||||
fn project(&self, p: Vec3, near: Option<[f64; 2]>, accept: f64) -> Option<[f64; 2]> {
|
||||
let mut q = self.surface.project(p, near, accept)?;
|
||||
if let Some(n) = near {
|
||||
for (d, x) in q.iter_mut().enumerate() {
|
||||
*x = uv::unwrap(*x, n[d], self.dom.period.get(d).copied().flatten());
|
||||
}
|
||||
}
|
||||
Some(q)
|
||||
}
|
||||
|
||||
fn sample(&self, curve: &CurveEval, t: f64, near: Option<[f64; 2]>, accept: 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 })
|
||||
}
|
||||
|
||||
/// 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);
|
||||
// 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);
|
||||
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)?;
|
||||
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.surface.pole(uv), 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 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<()> {
|
||||
let n = uses.len();
|
||||
let mut ends = Vec::with_capacity(n);
|
||||
for (c, fwd) in uses.iter() {
|
||||
let [a, b] = c.interval().ok_or_else(|| Error::Failed("pcurve interval".into()))?;
|
||||
let (Some(s), Some(e)) = (c.eval(a), c.eval(b)) else { return Err(Error::Failed("pcurve ends".into())) };
|
||||
ends.push(if *fwd { (s, e) } else { (e, s) });
|
||||
}
|
||||
for d in 0..2 {
|
||||
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 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];
|
||||
let mut free: Vec<usize> = Vec::new();
|
||||
let mut prev: Option<f64> = None;
|
||||
for (k, &i) in order.iter().enumerate() {
|
||||
let (Some(&(s, e)), Some(slot)) = (ends.get(i), shift.get_mut(i)) else { continue };
|
||||
if k > 0 && at_pole.get(i).copied().unwrap_or(false) {
|
||||
free.push(k);
|
||||
}
|
||||
*slot = prev.map_or(0.0, |p| ((p - s[d]) / period).round() * period);
|
||||
prev = Some(e[d] + *slot);
|
||||
}
|
||||
// After each later pole one period less, none or more: the right turn, then the
|
||||
// fewest moves.
|
||||
if !free.is_empty() && free.len() <= 6 {
|
||||
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);
|
||||
let mut extra = shift.clone();
|
||||
for (k, &i) in order.iter().enumerate() {
|
||||
if free.contains(&k) {
|
||||
let step = (code % 3) as f64 - 1.0;
|
||||
code /= 3;
|
||||
moved += usize::from(step != 0.0);
|
||||
add += step * period;
|
||||
}
|
||||
if let Some(slot) = extra.get_mut(i) {
|
||||
*slot += add;
|
||||
}
|
||||
}
|
||||
let right = area(uses, &extra, d) > 0.0;
|
||||
if best.as_ref().is_none_or(|(r, m, _)| (right && !r) || (right == *r && moved < *m)) {
|
||||
best = Some((right, moved, extra));
|
||||
}
|
||||
}
|
||||
if let Some((_, _, extra)) = best {
|
||||
shift = extra;
|
||||
}
|
||||
}
|
||||
let lo = self.dom.range.get(d).copied().flatten().map_or(0.0, |(a, b)| a.min(b));
|
||||
let low = uses.iter().zip(&shift).flat_map(|((c, _), s)| points(c).into_iter().map(move |q| q[d] + s)).fold(f64::INFINITY, f64::min);
|
||||
let base = if low.is_finite() { ((low - lo) / period + 1e-9).floor() * period } else { 0.0 };
|
||||
for (((c, _), s), (st, en)) in uses.iter_mut().zip(&shift).zip(ends.iter_mut()) {
|
||||
let by = s - base;
|
||||
if by != 0.0 {
|
||||
move_by(c, d, by);
|
||||
st[d] += by;
|
||||
en[d] += by;
|
||||
}
|
||||
}
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
}
|
||||
|
||||
/// 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).
|
||||
fn breaks(curve: &Curve, 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)| (k, (3 - degree.saturating_sub(m).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, t0, t1),
|
||||
},
|
||||
}
|
||||
}
|
||||
|
||||
/// 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();
|
||||
for ((c, fwd), s) in uses.iter().zip(shift) {
|
||||
let Some([a, b]) = c.interval() else { continue };
|
||||
for k in 0..16 {
|
||||
let f = k as f64 / 16.0;
|
||||
let t = if *fwd { a + (b - a) * f } else { b - (b - a) * f };
|
||||
if let Some(mut q) = c.eval(t) {
|
||||
if let Some(x) = q.get_mut(d) {
|
||||
*x += s;
|
||||
}
|
||||
ring.push(q);
|
||||
}
|
||||
}
|
||||
}
|
||||
let n = ring.len();
|
||||
(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 {
|
||||
Curve2::Spline { points, .. } => points.clone(),
|
||||
Curve2::RationalSpline { homogeneous, .. } => homogeneous.iter().map(|h| [h[0] / h[2], h[1] / h[2]]).collect(),
|
||||
}
|
||||
}
|
||||
|
||||
/// `c` moved by `by` along parameter `d`.
|
||||
fn move_by(c: &mut Curve2, d: usize, by: f64) {
|
||||
match c {
|
||||
Curve2::Spline { points, .. } => points.iter_mut().filter_map(|p| p.get_mut(d)).for_each(|x| *x += by),
|
||||
Curve2::RationalSpline { homogeneous, .. } => {
|
||||
for h in homogeneous {
|
||||
let w = h[2];
|
||||
if let Some(x) = h.get_mut(d) {
|
||||
*x += by * w;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
#[path = "pcurve_tests.rs"]
|
||||
mod tests;
|
||||
@@ -0,0 +1,310 @@
|
||||
use super::*;
|
||||
use crate::brep::Frame;
|
||||
use crate::brep::eval::{CurveEval, SurfaceEval};
|
||||
use std::f64::consts::{FRAC_PI_2, PI, TAU};
|
||||
|
||||
fn frame(origin: [f64; 3], x_axis: [f64; 3], y_axis: [f64; 3]) -> Frame {
|
||||
Frame { origin, x_axis, y_axis }
|
||||
}
|
||||
/// A frame about the axis through `origin` along +z, or −z when `down` (x stays +x).
|
||||
fn about_z(origin: [f64; 3], down: bool) -> Frame {
|
||||
frame(origin, [1., 0., 0.], if down { [0., -1., 0.] } else { [0., 1., 0.] })
|
||||
}
|
||||
fn one(curve: &Curve, interval: [f64; 2]) -> Vec<Vec<EdgeUse<'_>>> {
|
||||
vec![vec![EdgeUse { edge: 0, curve, interval, forward: true }]]
|
||||
}
|
||||
/// The largest distance between the lifted pcurve and the curve, at 65 points and each span's
|
||||
/// middle.
|
||||
fn lift_error(surface: &Surface, curve: &Curve, interval: [f64; 2], pc: &Curve2) -> f64 {
|
||||
let (s, c) = (SurfaceEval::new(surface).unwrap(), CurveEval::new(curve).unwrap());
|
||||
let [t0, t1] = interval;
|
||||
assert_eq!(pc.interval(), Some(interval), "on the edge's interval");
|
||||
let knots = match pc {
|
||||
Curve2::Spline { knots, .. } | Curve2::RationalSpline { knots, .. } => knots,
|
||||
};
|
||||
let mut ts: Vec<f64> = (0..=64).map(|k| t0 + (t1 - t0) * k as f64 / 64.).collect();
|
||||
ts.extend(knots.windows(2).filter(|k| k[0] < k[1]).map(|k| 0.5 * (k[0] + k[1])));
|
||||
ts.iter().map(|&t| s.point(pc.eval(t).unwrap()).unwrap().dist(c.point(t).unwrap())).fold(0., f64::max)
|
||||
}
|
||||
fn segment_points(pc: &Curve2) -> [[f64; 2]; 2] {
|
||||
match pc {
|
||||
Curve2::Spline { knots, points } if knots.len() == 4 && points.len() == 2 => [points[0], points[1]],
|
||||
other => panic!("not a degree-1 segment: {other:?}"),
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn lines_polylines_and_splines_on_planes_are_their_points_mapped() {
|
||||
let plane = Surface::PlaneThrough { origin: [1., 2., 3.], u_point: [3., 2., 3.], v_point: [1., 2., 7.] };
|
||||
let line = Curve::Line { start: [2., 2., 4.], end: [5., 2., 11.] };
|
||||
let pc = &pcurves(&plane, &one(&line, [0., 1.]), 1e-9).unwrap()[0][0];
|
||||
// Bit for bit: the ends lift to the line's own ends.
|
||||
let [a, b] = segment_points(pc);
|
||||
assert_eq!((a, b), ([0.5, 0.25], [2., 2.]));
|
||||
// On part of the line, extended past it.
|
||||
let pc = &pcurves(&plane, &one(&line, [-1., 3.]), 1e-9).unwrap()[0][0];
|
||||
assert!(lift_error(&plane, &line, [-1., 3.], pc) <= 1e-12);
|
||||
let poly = Curve::Polyline { points: vec![[1., 2., 3.], [3., 2., 3.], [3., 2., 7.], [5., 2., 9.]] };
|
||||
let pc = &pcurves(&plane, &one(&poly, [0.5, 3.]), 1e-9).unwrap()[0][0];
|
||||
let Curve2::Spline { knots, points } = pc else { panic!("a spline") };
|
||||
assert_eq!(knots, &[0.5, 0.5, 1., 2., 3., 3.]);
|
||||
assert_eq!(points, &[[0.5, 0.], [1., 0.], [1., 1.], [2., 1.5]]);
|
||||
assert!(lift_error(&plane, &poly, [0.5, 3.], pc) <= 1e-12);
|
||||
let w = std::f64::consts::FRAC_1_SQRT_2;
|
||||
let arc = Curve::RationalSpline {
|
||||
knots: vec![0., 0., 0., 1., 1., 1.],
|
||||
homogeneous: vec![[3., 2., 3., 1.], [3. * w, 2. * w, 5. * w, w], [1., 2., 5., 1.]],
|
||||
};
|
||||
let pc = &pcurves(&plane, &one(&arc, [0., 1.]), 1e-9).unwrap()[0][0];
|
||||
assert!(matches!(pc, Curve2::RationalSpline { .. }));
|
||||
assert!(lift_error(&plane, &arc, [0., 1.], pc) <= 1e-12);
|
||||
let cubic = Curve::Spline { knots: vec![0., 0., 0., 0., 2., 2., 2., 2.], points: vec![[1., 2., 3.], [4., 2., 3.], [2., 2., 6.], [5., 2., 8.]] };
|
||||
let pc = &pcurves(&plane, &one(&cubic, [0., 2.]), 1e-9).unwrap()[0][0];
|
||||
assert_eq!(pc, &Curve2::Spline { knots: vec![0., 0., 0., 0., 2., 2., 2., 2.], points: vec![[0., 0.], [1.5, 0.], [0.5, 0.75], [2., 1.25]] });
|
||||
// A piece of the spline is fitted instead (its knots would need cutting).
|
||||
let pc = &pcurves(&plane, &one(&cubic, [0.5, 1.5]), 1e-7).unwrap()[0][0];
|
||||
assert!(lift_error(&plane, &cubic, [0.5, 1.5], pc) <= 1e-7);
|
||||
}
|
||||
|
||||
/// A circle about the axis is a u-line: v fixed, u its angle plus a constant (minus where it
|
||||
/// turns the other way), on cylinders, cones, spheres and tori, in either direction.
|
||||
#[test]
|
||||
fn circles_about_the_axis_are_u_lines() {
|
||||
let o = [1., -2., 3.];
|
||||
let at = |h: f64| [1., -2., 3. + h];
|
||||
let tilt = 0.3f64;
|
||||
let rotated = |origin: [f64; 3], down: bool| {
|
||||
let (s, c) = tilt.sin_cos();
|
||||
frame(origin, [c, s, 0.], if down { [s, -c, 0.] } else { [-s, c, 0.] })
|
||||
};
|
||||
let half = 0.4f64;
|
||||
let cases = [
|
||||
(Surface::Cylinder { frame: about_z(o, false), radius: 5. }, 7., 5., 7.),
|
||||
(Surface::Cone { frame: about_z(o, false), radius_at_origin: 2., half_angle: half }, 4., 2. + 4. * half.tan(), 4.),
|
||||
(Surface::Sphere { frame: about_z(o, false), radius: 5. }, 3., 4., 3f64.atan2(4.)),
|
||||
(Surface::Sphere { frame: about_z(o, false), radius: 5. }, -3., 4., (-3f64).atan2(4.)),
|
||||
(Surface::Torus { frame: about_z(o, false), major_radius: 10., minor_radius: 2. }, 2f64.sqrt(), 10. - 2f64.sqrt(), 0.75 * PI),
|
||||
];
|
||||
for (surface, h, r, v) in cases {
|
||||
for down in [false, true] {
|
||||
let circle = Curve::Circle { frame: rotated(at(h), down), radius: r };
|
||||
for interval in [[0., TAU], [-1., 2.5], [4., 5.]] {
|
||||
let pc = &pcurves(&surface, &one(&circle, interval), 1e-9).unwrap()[0][0];
|
||||
let [a, b] = segment_points(pc);
|
||||
let sign = if down { -1. } else { 1. };
|
||||
// The u-line itself: the only freedom is the period the loop is moved by.
|
||||
let shift = a[0] - (tilt + sign * interval[0]);
|
||||
assert!((shift / TAU - (shift / TAU).round()).abs() < 1e-14, "{surface:?} {down}: {a:?}");
|
||||
assert!((b[0] - a[0] - sign * (interval[1] - interval[0])).abs() < 1e-13);
|
||||
assert!(a[1] == b[1] && (a[1] - v).abs() < 1e-14, "{surface:?}: v {a:?} {v}");
|
||||
let e = lift_error(&surface, &circle, interval, pc);
|
||||
assert!(e <= 1e-12, "{surface:?} {down} {interval:?}: {e}");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// A line along a ruling of a cylinder or a cone is a v-line: u fixed, v its height.
|
||||
#[test]
|
||||
fn rulings_are_v_lines() {
|
||||
let half = 0.3f64;
|
||||
let cylinder = Surface::Cylinder { frame: about_z([0.; 3], false), radius: 4. };
|
||||
let cone = Surface::Cone { frame: about_z([0., 0., 1.], false), radius_at_origin: 3., half_angle: half };
|
||||
let (s, c) = 2f64.sin_cos();
|
||||
let up = Curve::Line { start: [4. * c, 4. * s, -1.], end: [4. * c, 4. * s, 6.] };
|
||||
let down = Curve::Line { start: [4. * c, 4. * s, 6.], end: [4. * c, 4. * s, -1.] };
|
||||
let r = |z: f64| 3. + (z - 1.) * half.tan();
|
||||
let slant = Curve::Line { start: [r(2.) * c, r(2.) * s, 2.], end: [r(5.) * c, r(5.) * s, 5.] };
|
||||
for (surface, line, heights) in [(&cylinder, &up, [-1., 6.]), (&cylinder, &down, [6., -1.]), (&cone, &slant, [1., 4.])] {
|
||||
for interval in [[0., 1.], [0.25, 2.]] {
|
||||
let pc = &pcurves(surface, &one(line, interval), 1e-9).unwrap()[0][0];
|
||||
let [a, b] = segment_points(pc);
|
||||
assert_eq!(a[0], b[0]);
|
||||
assert!((a[0] - 2.).abs() < 1e-14);
|
||||
let height = |t: f64| heights[0] + (heights[1] - heights[0]) * t;
|
||||
assert!((a[1] - height(interval[0])).abs() < 1e-13 && (b[1] - height(interval[1])).abs() < 1e-13, "{a:?} {b:?}");
|
||||
assert!(lift_error(surface, line, interval, pc) <= 1e-12);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Curves with no closed form are fitted cubic B-splines within the tolerance: a tilted circle
|
||||
/// on a sphere, a plane's ellipse on a cylinder, a circle off the axis on a plane.
|
||||
#[test]
|
||||
fn other_curves_are_fitted_within_tolerance() {
|
||||
let (s, c) = 0.5f64.sin_cos();
|
||||
let sphere = Surface::Sphere { frame: about_z([0.; 3], false), radius: 6. };
|
||||
let great = Curve::Circle { frame: frame([0.; 3], [1., 0., 0.], [0., c, s]), radius: 6. };
|
||||
let cylinder = Surface::Cylinder { frame: about_z([0.; 3], false), radius: 2. };
|
||||
let ellipse = Curve::Ellipse {
|
||||
frame: frame([0., 0., 1.], [1. / (1f64 + 0.25).sqrt(), 0., 0.5 / (1f64 + 0.25).sqrt()], [0., 1., 0.]),
|
||||
radii: [2. * (1.25f64).sqrt(), 2.],
|
||||
};
|
||||
let plane = Surface::Plane { origin: [0.; 3], u_axis: [1., 0., 0.], v_axis: [0., 1., 0.] };
|
||||
let circle = Curve::Circle { frame: about_z([3., 1., 0.], false), radius: 2. };
|
||||
for (surface, curve, interval) in [(&sphere, &great, [0., TAU]), (&cylinder, &ellipse, [0.2, 6.]), (&plane, &circle, [0., 3.])] {
|
||||
for tol in [1e-3, 1e-6] {
|
||||
let pc = &pcurves(surface, &one(curve, interval), tol).unwrap()[0][0];
|
||||
let Curve2::Spline { knots, .. } = pc else { panic!("a B-spline") };
|
||||
assert_eq!(knots.len() - pc_count(pc) - 1, 3, "cubic");
|
||||
let e = lift_error(surface, curve, interval, pc);
|
||||
assert!(e <= tol, "{curve:?} at {tol}: {e}");
|
||||
}
|
||||
}
|
||||
}
|
||||
fn pc_count(pc: &Curve2) -> usize {
|
||||
match pc {
|
||||
Curve2::Spline { points, .. } => points.len(),
|
||||
Curve2::RationalSpline { homogeneous, .. } => homogeneous.len(),
|
||||
}
|
||||
}
|
||||
|
||||
/// The ends of each use in loop order and the loop's signed area in (u, v).
|
||||
fn walk(l: &[EdgeUse<'_>], pcs: &[Curve2]) -> (Vec<([f64; 2], [f64; 2])>, f64) {
|
||||
let mut ring = Vec::new();
|
||||
let mut ends = Vec::new();
|
||||
for (u, pc) in l.iter().zip(pcs) {
|
||||
let [a, b] = u.interval;
|
||||
let ts: Vec<f64> = (0..=16).map(|k| if u.forward { a + (b - a) * k as f64 / 16. } else { b - (b - a) * k as f64 / 16. }).collect();
|
||||
let pts: Vec<[f64; 2]> = ts.iter().map(|t| pc.eval(*t).unwrap()).collect();
|
||||
ends.push((pts[0], pts[16]));
|
||||
ring.extend_from_slice(&pts[..16]);
|
||||
}
|
||||
let n = ring.len();
|
||||
let area = (0..n).map(|i| ring[i][0] * ring[(i + 1) % n][1] - ring[(i + 1) % n][0] * ring[i][1]).sum::<f64>() * 0.5;
|
||||
(ends, area)
|
||||
}
|
||||
|
||||
/// A cylinder's wall between two full circles, with its seam line used twice: the seam's uses
|
||||
/// are one period apart, the loop closes and turns counter-clockwise, all in the base period.
|
||||
#[test]
|
||||
fn a_seam_used_twice_gets_pcurves_a_period_apart() {
|
||||
let cylinder = Surface::Cylinder { frame: about_z([0.; 3], false), radius: 3. };
|
||||
let seam = Curve::Line { start: [3., 0., 0.], end: [3., 0., 4.] };
|
||||
let low = Curve::Circle { frame: about_z([0.; 3], false), radius: 3. };
|
||||
let high = Curve::Circle { frame: about_z([0., 0., 4.], false), radius: 3. };
|
||||
let full = [0., TAU];
|
||||
let l = vec![
|
||||
EdgeUse { edge: 0, curve: &low, interval: full, forward: true },
|
||||
EdgeUse { edge: 1, curve: &seam, interval: [0., 1.], forward: true },
|
||||
EdgeUse { edge: 2, curve: &high, interval: full, forward: false },
|
||||
EdgeUse { edge: 1, curve: &seam, interval: [0., 1.], forward: false },
|
||||
];
|
||||
let pcs = &pcurves(&cylinder, 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 > 0., "counter-clockwise: {area}");
|
||||
let (a, b) = (segment_points(&pcs[1]), segment_points(&pcs[3]));
|
||||
assert!(((a[0][0] - b[0][0]).abs() - TAU).abs() < 1e-12 && a[0][1] == b[0][1]);
|
||||
for pc in pcs {
|
||||
for p in segment_points(pc) {
|
||||
assert!((-1e-12..=TAU + 1e-12).contains(&p[0]), "base period: {p:?}");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// A sphere bounded by one meridian used twice from pole to pole: each use runs straight into
|
||||
/// the poles, the two a period apart, and the loop runs along the poles' lines.
|
||||
#[test]
|
||||
fn meridians_run_straight_into_the_poles() {
|
||||
let sphere = Surface::Sphere { frame: about_z([0.; 3], false), radius: 6. };
|
||||
let meridian = Curve::Circle { frame: frame([0.; 3], [1., 0., 0.], [0., 0., 1.]), radius: 6. };
|
||||
let half = [-FRAC_PI_2, FRAC_PI_2];
|
||||
let l = vec![
|
||||
EdgeUse { edge: 0, curve: &meridian, interval: half, forward: false },
|
||||
EdgeUse { edge: 0, curve: &meridian, interval: half, forward: true },
|
||||
];
|
||||
for tol in [1e-3, 1e-6] {
|
||||
let pcs = &pcurves(&sphere, std::slice::from_ref(&l), tol).unwrap()[0];
|
||||
for pc in pcs {
|
||||
assert!(lift_error(&sphere, &meridian, half, pc) <= tol);
|
||||
for k in 0..=16 {
|
||||
let q = pc.eval(-FRAC_PI_2 + PI * k as f64 / 16.).unwrap();
|
||||
let u = pcs[0].eval(0.).unwrap()[0];
|
||||
assert!((q[0] - u).abs() < 1e-9 || ((q[0] - u).abs() - TAU).abs() < 1e-9, "straight: {q:?}");
|
||||
}
|
||||
}
|
||||
let (a, b) = (pcs[0].eval(0.).unwrap(), pcs[1].eval(0.).unwrap());
|
||||
assert!(((a[0] - b[0]).abs() - TAU).abs() < 1e-9 && (a[1] - b[1]).abs() < 1e-12, "{a:?} {b:?}");
|
||||
let (_, area) = walk(&l, pcs);
|
||||
assert!(area > 0., "counter-clockwise: {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]
|
||||
fn intersection_curves_lie_on_both_surfaces() {
|
||||
// The cylinder x² + y² = 25 cut by the plane z = 20 + x tan 0.4: an ellipse.
|
||||
let cylinder = Surface::Cylinder { frame: about_z([0.; 3], false), radius: 5. };
|
||||
let k = 0.4f64.tan();
|
||||
let plane = Surface::Plane { origin: [0., 0., 20.], u_axis: [1., 0., k], v_axis: [0., 1., 0.] };
|
||||
// A coarse leader through points of the curve, its chords up to 0.03 mm off it.
|
||||
let leader: Vec<[f64; 3]> = (0..=12)
|
||||
.map(|i| {
|
||||
let a = 0.2 + 2.5 * i as f64 / 12.;
|
||||
[5. * a.cos(), 5. * a.sin(), 20. + 5. * a.cos() * k]
|
||||
})
|
||||
.collect();
|
||||
let curve = Curve::SurfaceIntersection {
|
||||
surfaces: [Box::new(cylinder.clone()), Box::new(plane.clone())],
|
||||
leader: Box::new(Curve::Polyline { points: leader }),
|
||||
tolerance: 0.2,
|
||||
};
|
||||
let eval = CurveEval::new(&curve).unwrap();
|
||||
for i in 0..=48 {
|
||||
let p = eval.point(12. * i as f64 / 48.).unwrap();
|
||||
assert!((p.x.hypot(p.y) - 5.).abs() < 1e-11 && (p.z - 20. - p.x * k).abs() < 1e-11, "{i}: {p:?}");
|
||||
}
|
||||
for surface in [&cylinder, &plane] {
|
||||
let pc = &pcurves(surface, &one(&curve, [0., 12.]), 1e-6).unwrap()[0][0];
|
||||
assert!(lift_error(surface, &curve, [0., 12.], pc) <= 1e-6);
|
||||
}
|
||||
}
|
||||
|
||||
/// Hostile records give errors, never a panic.
|
||||
#[test]
|
||||
fn hostile_records_give_errors() {
|
||||
let plane = Surface::Plane { origin: [0.; 3], u_axis: [1., 0., 0.], v_axis: [0., 1., 0.] };
|
||||
let line = Curve::Line { start: [0.; 3], end: [1., 0., 0.] };
|
||||
for tol in [0., -1., f64::NAN, f64::INFINITY, 1e9] {
|
||||
assert!(matches!(pcurves(&plane, &one(&line, [0., 1.]), tol), Err(Error::Invalid(_))), "{tol}");
|
||||
}
|
||||
for interval in [[1., 1.], [2., 1.], [f64::NAN, 1.], [0., f64::INFINITY], [0., 1e13]] {
|
||||
assert!(matches!(pcurves(&plane, &one(&line, interval), 1e-6), Err(Error::Invalid(_))), "{interval:?}");
|
||||
}
|
||||
let nan = Curve::Line { start: [f64::NAN, 0., 0.], end: [1., 0., 0.] };
|
||||
assert!(matches!(pcurves(&plane, &one(&nan, [0., 1.]), 1e-6), Err(Error::Invalid(_))));
|
||||
let skew = Surface::Cylinder { frame: frame([0.; 3], [1., 0., 0.], [0.5, 1., 0.]), radius: 1. };
|
||||
assert!(matches!(pcurves(&skew, &one(&line, [0., 1.]), 1e-6), Err(Error::Invalid(_))));
|
||||
let flat = Surface::Plane { origin: [0.; 3], u_axis: [1., 0., 0.], v_axis: [2., 0., 0.] };
|
||||
assert!(matches!(pcurves(&flat, &one(&line, [0., 1.]), 1e-6), Err(Error::Invalid(_))));
|
||||
let collapsed = Curve::Circle { frame: about_z([0.; 3], false), radius: 0. };
|
||||
assert!(matches!(pcurves(&plane, &one(&collapsed, [0., 1.]), 1e-6), Err(Error::Invalid(_))));
|
||||
// Too many uses, refused before any work.
|
||||
let many = vec![vec![EdgeUse { edge: 0, curve: &line, interval: [0., 1.], forward: true }; MAX_USES + 1]];
|
||||
assert!(matches!(pcurves(&plane, &many, 1e-6), Err(Error::TooLarge(_))));
|
||||
// A polyline of many points on a curved surface: too many spans for one pcurve.
|
||||
let cylinder = Surface::Cylinder { frame: about_z([0.; 3], false), radius: 1. };
|
||||
let n = 5000;
|
||||
let zigzag = Curve::Polyline { points: (0..n).map(|i| [(i as f64 * 0.01).cos(), (i as f64 * 0.01).sin(), (i % 2) as f64]).collect() };
|
||||
assert!(matches!(pcurves(&cylinder, &one(&zigzag, [0., (n - 1) as f64]), 1e-6), Err(Error::TooLarge(_))));
|
||||
// A degree the evaluators do not take.
|
||||
let high = Curve::Spline { knots: [vec![0.; 10], vec![1.; 10]].concat(), points: (0..10).map(|i| [i as f64, 0., 0.]).collect() };
|
||||
assert!(matches!(pcurves(&cylinder, &one(&high, [0., 1.]), 1e-6), Err(Error::Failed(_))));
|
||||
// Off the surface, and a curve on two tangent surfaces (no single point).
|
||||
assert!(matches!(pcurves(&cylinder, &one(&line, [0., 1.]), 1e-6), Err(Error::Failed(_))));
|
||||
let touching = Curve::SurfaceIntersection {
|
||||
surfaces: [Box::new(cylinder.clone()), Box::new(Surface::Plane { origin: [1., 0., 0.], u_axis: [0., 1., 0.], v_axis: [0., 0., 1.] })],
|
||||
leader: Box::new(Curve::Polyline { points: vec![[1., 0., 0.], [1., 0., 1.]] }),
|
||||
tolerance: 1e-3,
|
||||
};
|
||||
assert!(pcurves(&cylinder, &one(&touching, [0., 1.]), 1e-6).is_err());
|
||||
// Extreme but valid: a torus as wide as the schema allows and thin as it allows.
|
||||
let torus = Surface::Torus { frame: about_z([0.; 3], false), major_radius: 1e12, minor_radius: 1e-12 };
|
||||
let far = Curve::Circle { frame: about_z([0.; 3], false), radius: 1e12 };
|
||||
let _ = pcurves(&torus, &one(&far, [0., 1.]), 1e-6);
|
||||
}
|
||||
@@ -6,7 +6,7 @@ use super::EvalSurface;
|
||||
use crate::Vec3;
|
||||
|
||||
/// The parameter box of a surface and its periods.
|
||||
#[derive(Clone, Copy, Debug)]
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub struct Domain {
|
||||
/// Bounded ranges (`None`: unbounded, as a plane's).
|
||||
pub range: [Option<(f64, f64)>; 2],
|
||||
|
||||
Reference in New Issue
Block a user