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

Drawing W06a: exact orthographic hidden lines for the analytic subset (#2455) (#2457)

Browse Source 
* Drawing W06a: exact orthographic hidden lines for the analytic subset (#2455)

crates/geom/src/hlr (L0): prepare once per body, then hidden_lines for a
scene of placed bodies. Edges and face silhouettes are projected from their
exact curves, split where projections cross, and each piece is classified by
rays toward the viewer against the trimmed faces: visible, hidden or
explicitly uncertain, never guessed. Planes, cylinders, cones, spheres and
tori are supported whether stored as analytic records, revolutions or spline
patches; freeform faces are reported and make what they could hide
uncertain. Screen geometry is Seg2 with the W01a accuracy vocabulary: exact
lines and conic arcs where the projection is exact, bounded chords
otherwise. docs/drawing-hlr.md holds the contract and its limits.

* geom/hlr: handle the brep records main gained in #3093

Main's exact model storage (#3093) added Curve::SurfaceIntersection,
Surface::PlaneThrough and Surface::PlacedRevolution, and the HLR matches on
the records no longer compiled. A plane through three points is a plane: it
is trimmed and intersected as one, and a box stored that way draws the same
bits as one with axis-stored planes. The procedural intersection curve and
the placed revolution (whose affine placement may scale it unevenly) are
refused by prepare as not supported yet, until the kernel adapter handles
them. The records' positions move with the body's centre like the others.

* geom/hlr: pad plane faces' boxes, so a plane record off by rounding still hides

Plane faces took the box of their boundary chords unpadded, flat along the
normal. When the plane record and its edges disagree by rounding (a unit cube's
top plane one ulp above its edges, a 100 mm cube's 1e-11), the ray toward the
viewer reached the plane just outside that box, the slab test culled the hit,
and the back edges came out Visible with no doubt. Every face's box is now
padded by 1e-6 of its size, as curved faces' boxes already were. Reported by
DevonHope on #2457.

* geom/hlr: limits from review: trim samples, scene extent, crossing pairs, surface groups

From DevonHope's review of #2457:
- Options::trim_samples must be 1 to MAX_TRIM_SAMPLES (4096).
- Placed bodies must stay within MAX_COORDINATE (1e12 mm, as section and brep
  bound stored numbers). Placements 1e200 apart made the scene's size infinite
  and every tolerance taken from it with it, so everything came out Visible.
- The crossing search counts its candidate pairs before it starts and meets
  each pair only in the first grid cell the two segments share, instead of
  keeping a set of the pairs met, which could outgrow wasm's 4 GB before the
  500 million pair budget tripped. The same pairs in the same order.
- prepare groups faces on one quadric or torus by comparing each face with
  the faces of its kind whose size agrees to the tolerance (sorted), not with
  every face; at most 50 million comparisons. prepare_cancellable polls a
  cancel callback between edges, faces and comparisons.

* geom/hlr: boundary chords on a grid, and no fine pass over coincident curved chords

The cost items of DevonHope's review of #2457:
- Every ray hit asked how far it lies from its face's boundary by scanning all
  of the face's chords (4096 per circle; a plate's face with 256 holes has
  about a million). The chords are now bucketed on a uniform grid over the
  boundary's box, cells at least as long as the longest chord, and the test
  looks only at the cells around the point. It finds exactly the chords the
  scan finds.
- Chords of two curved curves that lie on each other end to end (aligned holes
  seen along their axes) were cut into 16 by 16 pieces and intersected, to
  find crossings that are only noise within the curves' tolerance. They are
  now skipped; visibility changes where a coincident run starts or ends, and
  chords that only partly overlap still find those.
- Fixtures: a plate with holes, coaxial cylinders seen along their axis, and
  300 placements.

* geom/hlr: pad boxes by the body's coordinates too; coincident chords only to rounding

Review of the previous commits found two holes:
- The padding grew with the face's size, but rounding grows with where the body
  was modelled: a 0.001 mm cube 1e7 mm out, its top plane one step (1.9e-9 mm)
  above its edges, still had its hit culled. The pad now adds 1e-12 of the
  body's offset.
- Two chords were taken for coincident when their ends lay within each other's
  sag, which in a coarse scene can hide real crossings: an arc can bulge across
  a line 1 mm off its chord twice. Only chords that coincide to rounding are
  skipped now; identical projections, the costly case, still are.

* geom/hlr: coincident curves meeting at shared sample points do not cross

Timing the review's scenes showed that skipping coincident chords was not
enough: chords of two coincident curves also touch end to end at every sample
point they share, and each touch was a crossing to refine and split at. 64
coaxial cylinders seen along their axis gave 8.3 million crossings (19 s), 64
cylinders one behind another along the line of sight 4.2 million (214 s;
before these commits the crossing budget refused it). A touch at a shared
sample point is now dropped when the two curves also coincide, to rounding, on
both sides of it. Where a coincident run starts or ends one side differs, and
the touch is kept.

* geom/hlr: the own-face margin allows for rounding where a body was modelled

The far 0.001 mm cube of the previous commit showed a second effect of records
rounded at large coordinates: its top plane, one ulp (1.9e-9 mm) above its
edges, hid the top edges themselves, because a ray leaving an edge met its own
face 3.2e-9 mm along, past the own-face margin of 1e-6 of the face (1.4e-9).
The margin now also covers 1e-14 of the body's offset (about 45 ulps of its
coordinates), which only matters for parts far smaller than their distance
from the origin.

* geom/hlr: no parity test for a hit on the boundary, where it is uncertain anyway

* docs/drawing-hlr: the fixtures of the review's second round

---------

Co-authored-by: MateuszSychXpl <151755914+MateuszSychXpl@users.noreply.github.com>
This commit is contained in:
Mateusz Sych
2026-10-27 05:20:17 +01:00
committed by GitHub
parent 3575c99750
commit ebafe07487
8 changed files with 6725 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
crates/geom/src/hlr/curve.rs
+371
View File
@@ -0,0 +1,371 @@
//! Evaluation of portable B-rep curves: lines, circle and ellipse arcs, polynomial and rational
//! splines (fitted intersection records are ordinary splines with a stated fitting tolerance).
use super::HlrError;
use crate::Vec3;
use crate::brep::{Curve, Frame};
use crate::section::CurveAccuracy;
#[derive(Clone, Debug)]
pub(crate) enum C3 {
Line {
a: Vec3,
b: Vec3,
},
/// origin + x*rx*cos(t) + y*ry*sin(t), t in `range`.
Arc {
o: Vec3,
x: Vec3,
y: Vec3,
r: [f64; 2],
range: [f64; 2],
},
/// Homogeneous control points `[x*w, y*w, z*w, w]` (w = 1 for polynomial splines), and the
/// same curve as Bézier pieces (each piece's parameter span and control points; empty when
/// it could not be split).
Spline {
degree: usize,
knots: Vec<f64>,
points: Vec<[f64; 4]>,
pieces: Vec<([f64; 2], Vec<[f64; 4]>)>,
},
Polyline(Vec<Vec3>),
}
pub(crate) fn v3(p: [f64; 3]) -> Vec3 {
Vec3::new(p[0], p[1], p[2])
}
pub(crate) fn frame_axes(f: &Frame) -> (Vec3, Vec3, Vec3) {
let (x, y) = (v3(f.x_axis), v3(f.y_axis));
(v3(f.origin), x, y)
}
impl C3 {
pub fn new(c: &Curve) -> Result<(C3, CurveAccuracy), HlrError> {
Ok(match c {
Curve::Line { start, end } => (C3::Line { a: v3(*start), b: v3(*end) }, CurveAccuracy::Analytic),
Curve::CircleArc { frame, radius, angle } => {
let (o, x, y) = frame_axes(frame);
(C3::Arc { o, x, y, r: [*radius, *radius], range: *angle }, CurveAccuracy::Analytic)
}
Curve::EllipseArc { frame, radii, angle } => {
let (o, x, y) = frame_axes(frame);
(C3::Arc { o, x, y, r: *radii, range: *angle }, CurveAccuracy::Analytic)
}
Curve::Polyline { points } => (C3::Polyline(points.iter().copied().map(v3).collect()), CurveAccuracy::Analytic),
Curve::Spline { knots, points } => (spline(knots, points.iter().map(|p| [p[0], p[1], p[2], 1.]).collect())?, CurveAccuracy::Analytic),
Curve::RationalSpline { knots, homogeneous } => (spline(knots, homogeneous.clone())?, CurveAccuracy::Analytic),
Curve::Intersection { knots, points, fitting_tolerance, .. } => (
spline(knots, points.iter().map(|p| [p[0], p[1], p[2], 1.]).collect())?,
CurveAccuracy::Bounded { max_deviation_mm: *fitting_tolerance },
),
// Procedural: the curve is traced on demand from its leader, which carries no bound.
Curve::SurfaceIntersection { .. } => return Err(HlrError::Unsupported("surface intersection curve")),
})
}
pub fn domain(&self) -> [f64; 2] {
match self {
C3::Line { .. } => [0., 1.],
C3::Arc { range, .. } => *range,
C3::Spline { degree, knots, points, .. } => [knots.get(*degree).copied().unwrap_or(0.), knots.get(points.len()).copied().unwrap_or(1.)],
C3::Polyline(p) => [0., p.len().saturating_sub(1) as f64],
}
}
pub fn at(&self, t: f64) -> Vec3 {
match self {
C3::Line { a, b } => *a + (*b - *a) * t,
C3::Arc { o, x, y, r, .. } => *o + *x * (r[0] * t.cos()) + *y * (r[1] * t.sin()),
C3::Spline { degree, knots, points, .. } => de_boor(*degree, knots, points, t),
C3::Polyline(p) => {
let n = p.len();
if n == 0 {
return Vec3::ZERO;
}
let t = t.clamp(0., (n - 1) as f64);
let i = (t.floor() as usize).min(n.saturating_sub(2));
match (p.get(i), p.get(i + 1)) {
(Some(a), Some(b)) => a.lerp(*b, t - i as f64),
(Some(a), None) => *a,
_ => Vec3::ZERO,
}
}
}
}
/// First derivative by a central difference scaled to the domain.
pub fn d1(&self, t: f64) -> Vec3 {
match self {
C3::Line { a, b } => *b - *a,
C3::Arc { x, y, r, .. } => *x * (-r[0] * t.sin()) + *y * (r[1] * t.cos()),
_ => {
let [t0, t1] = self.domain();
let h = (t1 - t0) * 1e-6;
let (a, b) = ((t - h).max(t0), (t + h).min(t1));
if b - a <= 0. { Vec3::ZERO } else { (self.at(b) - self.at(a)) / (b - a) }
}
}
}
/// Parameters that split the curve into pieces each turning at most `max_turn` radians or
/// spanning at most a quarter of the domain, plus all interior knots.
pub fn breaks(&self) -> Vec<f64> {
let [t0, t1] = self.domain();
let mut out = vec![t0, t1];
match self {
C3::Spline { knots, .. } => out.extend(knots.iter().copied().filter(|k| *k > t0 && *k < t1)),
C3::Polyline(p) => out.extend((1..p.len().saturating_sub(1)).map(|i| i as f64)),
C3::Arc { .. } => {
let n = ((t1 - t0) / (std::f64::consts::PI / 8.)).ceil().max(1.) as usize;
out.extend((1..n).map(|i| t0 + (t1 - t0) * i as f64 / n as f64));
}
C3::Line { .. } => {}
}
out.sort_by(f64::total_cmp);
// Only equal knots merge: a span narrower than rounding can still carry a whole arc.
out.dedup();
out
}
}
fn spline(knots: &[f64], points: Vec<[f64; 4]>) -> Result<C3, HlrError> {
let degree = knots.len().checked_sub(points.len() + 1).ok_or(HlrError::Unsupported("spline knot count"))?;
if degree == 0 || degree > 16 || points.is_empty() {
return Err(HlrError::Unsupported("spline degree"));
}
let pieces = bezier_pieces(degree, knots, &points);
Ok(C3::Spline { degree, knots: knots.to_vec(), points, pieces })
}
/// A spline split into Bézier pieces: clamped, then every inner knot raised to the degree
/// (NURBS Book §5.2); empty when that fails.
fn bezier_pieces(p: usize, knots: &[f64], points: &[[f64; 4]]) -> Vec<([f64; 2], Vec<[f64; 4]>)> {
use crate::bspline::{HPoint, clamp_row, insert_knot};
let row: Vec<HPoint> = points.iter().map(|h| HPoint { x: h[0], y: h[1], z: h[2], w: h[3] }).collect();
let Ok((mut k, mut c)) = clamp_row(p, knots, &row) else { return Vec::new() };
let (Some(&lo), Some(&hi)) = (k.get(p), k.get(c.len())) else { return Vec::new() };
let mut inner: Vec<f64> = k.iter().copied().filter(|x| *x > lo && *x < hi).collect();
inner.dedup();
for u in &inner {
let mut guard = 0;
while k.iter().filter(|x| *x == u).count() < p {
guard += 1;
let Ok((nk, nc)) = insert_knot(p, &k, &c, *u) else { return Vec::new() };
if guard > p + 1 {
return Vec::new();
}
(k, c) = (nk, nc);
}
}
let breaks: Vec<f64> = std::iter::once(lo).chain(inner).chain(std::iter::once(hi)).collect();
let mut out = Vec::with_capacity(breaks.len());
for (j, w) in breaks.windows(2).enumerate() {
let (Some(&s0), Some(&s1)) = (w.first(), w.get(1)) else { return Vec::new() };
let Some(ctrl) = c.get(j * p..=j * p + p) else { return Vec::new() };
out.push(([s0, s1], ctrl.iter().map(|h| [h.x, h.y, h.z, h.w]).collect()));
}
out
}
/// The Bernstein coefficients of the product of two polynomials given by theirs.
fn bernstein_product(a: &[f64], b: &[f64]) -> Vec<f64> {
use crate::bspline::binomial;
let (m, n) = (a.len().saturating_sub(1), b.len().saturating_sub(1));
(0..=m + n)
.map(|k| {
let terms = (k.saturating_sub(n)..=k.min(m))
.map(|i| binomial(m, i) * binomial(n, k - i) * a.get(i).copied().unwrap_or(0.) * b.get(k - i).copied().unwrap_or(0.));
terms.sum::<f64>() / binomial(m + n, k)
})
.collect()
}
/// The sign a polynomial in Bernstein form keeps over its interval: +1, -1, 0 when it is zero
/// (to rounding) all along, `None` when that cannot be certified. All coefficients of one sign
/// mean the polynomial has it (their hull holds its graph); otherwise the interval is halved.
fn definite_sign(c: &[f64], depth: u32) -> Option<f64> {
let scale = c.iter().fold(0f64, |m, x| m.max(x.abs()));
let eps = 1e-12 * scale;
if c.iter().all(|x| x.abs() <= eps) {
return Some(0.);
}
if c.iter().all(|x| *x >= -eps) {
return Some(1.);
}
if c.iter().all(|x| *x <= eps) {
return Some(-1.);
}
if depth >= 16 {
return None;
}
// de Casteljau at the middle.
let (mut left, mut right, mut pts) = (Vec::with_capacity(c.len()), Vec::with_capacity(c.len()), c.to_vec());
while let (Some(&f), Some(&l)) = (pts.first(), pts.last()) {
left.push(f);
right.push(l);
pts = pts.windows(2).map(|w| 0.5 * (w.first().copied().unwrap_or(0.) + w.get(1).copied().unwrap_or(0.))).collect();
}
right.reverse();
match (definite_sign(&left, depth + 1)?, definite_sign(&right, depth + 1)?) {
(a, b) if a == b || b == 0. => Some(a),
(0., b) => Some(b),
_ => None,
}
}
impl C3 {
/// The way a spline's part `r` turns about the centre `o` in the plane of `u` and `v`:
/// +1 or -1 when its angular speed keeps that sign all along (certified on each Bézier piece
/// from the Bernstein coefficients of `X Y' - Y X'`, which has the speed's sign as the
/// weights are positive), 0 when it does not turn, `None` when it turns back or this cannot
/// be certified.
pub fn turn_sign(&self, r: [f64; 2], o: Vec3, u: Vec3, v: Vec3) -> Option<f64> {
use crate::bspline::{HPoint, bezier_segment};
let C3::Spline { pieces, .. } = self else { return None };
let mut sign = 0f64;
for (span, ctrl) in pieces {
let (a, b) = (r[0].max(span[0]), r[1].min(span[1]));
let len = span[1] - span[0];
if !(b > a && len > 0.) {
continue;
}
let (fa, fb) = (((a - span[0]) / len).clamp(0., 1.), ((b - span[0]) / len).clamp(0., 1.));
if !(fb > fa) {
continue;
}
let h: Vec<HPoint> = ctrl.iter().map(|q| HPoint { x: q[0], y: q[1], z: q[2], w: q[3] }).collect();
let seg = bezier_segment(&h, fa, fb).ok()?;
if seg.iter().any(|q| !(q.w > 0.)) {
return None;
}
let rel = |q: &HPoint| Vec3::new(q.x, q.y, q.z) - o * q.w;
let x: Vec<f64> = seg.iter().map(|q| rel(q).dot(u)).collect();
let y: Vec<f64> = seg.iter().map(|q| rel(q).dot(v)).collect();
let p = x.len().saturating_sub(1);
if p == 0 {
return None;
}
let d = |c: &[f64]| {
c.windows(2).map(|w| p as f64 * (w.get(1).copied().unwrap_or(0.) - w.first().copied().unwrap_or(0.))).collect::<Vec<f64>>()
};
let n: Vec<f64> = bernstein_product(&x, &d(&y)).iter().zip(bernstein_product(&y, &d(&x))).map(|(a, b)| a - b).collect();
let s = definite_sign(&n, 0)?;
if s != 0. {
if sign != 0. && s != sign {
return None;
}
sign = s;
}
}
Some(sign)
}
/// Points whose hull holds the spline's part `r`: the projected control points of its
/// Bézier pieces cut to `r` (de Casteljau). They close in on the curve as `r` shrinks.
/// `None` for other curves, or when a weight is not positive.
pub fn hull(&self, r: [f64; 2]) -> Option<Vec<Vec3>> {
use crate::bspline::{HPoint, bezier_segment};
let C3::Spline { pieces, .. } = self else { return None };
let mut out = Vec::new();
for (span, ctrl) in pieces {
let (a, b) = (r[0].max(span[0]), r[1].min(span[1]));
let len = span[1] - span[0];
if !(b > a && len > 0.) {
continue;
}
let (fa, fb) = (((a - span[0]) / len).clamp(0., 1.), ((b - span[0]) / len).clamp(0., 1.));
if !(fb > fa) {
continue;
}
let h: Vec<HPoint> = ctrl.iter().map(|q| HPoint { x: q[0], y: q[1], z: q[2], w: q[3] }).collect();
for q in bezier_segment(&h, fa, fb).ok()? {
if !(q.w > 0.) {
return None;
}
out.push(q.project()?);
}
}
(!out.is_empty()).then_some(out)
}
}
/// De Boor evaluation of a (rational) B-spline at `t`, clamped to the active domain.
pub(crate) fn de_boor(p: usize, knots: &[f64], ctrl: &[[f64; 4]], t: f64) -> Vec3 {
let n = ctrl.len();
let (Some(&lo), Some(&hi)) = (knots.get(p), knots.get(n)) else { return Vec3::ZERO };
let t = t.clamp(lo, hi);
// Span k with knots[k] <= t < knots[k+1], last nonempty span at the end.
let mut k = p;
while k + 1 < n && knots.get(k + 1).is_some_and(|u| *u <= t) {
k += 1;
}
let mut d: Vec<[f64; 4]> = (0..=p).map(|j| ctrl.get(j + k - p).copied().unwrap_or([0.; 4])).collect();
for r in 1..=p {
for j in (r..=p).rev() {
let (a, b) = (knots.get(j + k - p).copied().unwrap_or(0.), knots.get(j + 1 + k - r).copied().unwrap_or(0.));
let alpha = if b - a == 0. { 0. } else { (t - a) / (b - a) };
let (prev, cur) = (d.get(j - 1).copied().unwrap_or([0.; 4]), d.get(j).copied().unwrap_or([0.; 4]));
if let Some(slot) = d.get_mut(j) {
for c in 0..4 {
slot[c] = (1. - alpha) * prev[c] + alpha * cur[c];
}
}
}
}
let h = d.get(p).copied().unwrap_or([0.; 4]);
if h[3] == 0. { Vec3::ZERO } else { Vec3::new(h[0] / h[3], h[1] / h[3], h[2] / h[3]) }
}
/// The circle a curve lies on, if it is one: centre, unit normal, radius. Checked at many
/// parameters against `tol` (relative to the radius), so rational circle splines are recognised.
pub(crate) fn as_circle(c: &C3, tol: f64) -> Option<(Vec3, Vec3, f64)> {
match c {
C3::Arc { o, x, y, r, .. } if (r[0] - r[1]).abs() <= tol * r[0] => Some((*o, x.cross(*y).normalized()?, r[0])),
C3::Spline { degree, points, .. } => {
let [t0, t1] = c.domain();
let at = |f: f64| c.at(t0 + (t1 - t0) * f);
let (a, b, m) = (at(0.), at(1.), at(0.5));
// Circumcentre of a, m and a third point (b may equal a for a full circle).
let third = if a.dist(b) < 1e-9 * (1. + a.len()) { at(0.25) } else { b };
let (u, v) = (m - a, third - a);
let n = u.cross(v);
let n2 = n.len2();
if n2 <= 1e-24 {
return None;
}
let centre = a + (v * u.len2() - u * v.len2()).cross(n) * (1. / (2. * n2));
let normal = n.normalized()?;
let r = centre.dist(a);
// On one knot span of a (rational) spline of degree p, f = w^2 (|P - c|^2 - r^2) is
// a polynomial of degree 2p. At 2p + 3 Chebyshev points of the span it is at most
// w_max^2 2 r d (d the largest distance from the circle there), so everywhere at most
// 4 times that (their Lebesgue constant is below 4), and as w >= w_min the distance
// anywhere is at most 4 (w_max / w_min)^2 d: weights near zero cannot hide a swerve.
let weights = points.iter().map(|h| h[3]);
let (w_min, w_max) = weights.fold((f64::MAX, f64::MIN), |(lo, hi), w| (lo.min(w), hi.max(w)));
if !(w_min > 0. && w_max.is_finite()) {
return None;
}
let off = |p: Vec3| (p.dist(centre) - r).abs().max((p - centre).dot(normal).abs());
let n = (2 * degree + 3).max(9);
let mut d = 0f64;
for w in c.breaks().windows(2) {
let (Some(&s0), Some(&s1)) = (w.first(), w.get(1)) else { continue };
let (mid, half) = (0.5 * (s0 + s1), 0.5 * (s1 - s0));
d = d.max(off(c.at(s0))).max(off(c.at(s1)));
for k in 0..n {
d = d.max(off(c.at(mid + half * (std::f64::consts::PI * (2 * k + 1) as f64 / (2 * n) as f64).cos())));
}
}
// Relative to the curve's own size when that is smaller than the radius: a nearly
// straight curve fits a huge circle, whose radius says nothing about the error.
let pts: Vec<Vec3> = points.iter().filter(|h| h[3] != 0.).map(|h| Vec3::new(h[0] / h[3], h[1] / h[3], h[2] / h[3])).collect();
let size = pts.iter().skip(1).fold(pts.first().map(|p| (*p, *p)), |acc, p| acc.map(|(lo, hi)| (lo.min(*p), hi.max(*p))));
let size = size.map_or(r, |(lo, hi)| (hi - lo).len());
(4. * (w_max / w_min).powi(2) * d <= tol * r.min(size)).then_some((centre, normal, r))
}
_ => None,
}
}
crates/geom/src/hlr/mod.rs
+1915
View File
File diff suppressed because it is too large Load Diff
crates/geom/src/hlr/silhouette.rs
+347
View File
@@ -0,0 +1,347 @@
//! Silhouettes of curved faces: the curves where the surface normal is perpendicular to the
//! direction of sight. Recognised cylinders and cones give exact lines and spheres an exact
//! circle; other surfaces are contoured on a chart grid, each point refined onto the curve,
//! and the deviation of the chords from the curve is measured.
use super::curve::C3;
use super::surface::{Analytic, S3};
use super::trim::Trim;
use crate::Vec3;
use crate::section::CurveAccuracy;
const TAU: f64 = std::f64::consts::TAU;
pub(crate) struct Sil {
pub curve: C3,
pub acc: CurveAccuracy,
/// Parameter ranges of `curve` that lie inside the face.
pub ranges: Vec<[f64; 2]>,
}
fn perp(a: Vec3) -> Vec3 {
a.any_perp().normalized().unwrap_or(Vec3::new(1., 0., 0.))
}
/// Axial extent of the faces' boundaries relative to `p` along `a`.
fn extent(parts: &[(&S3, &Trim)], p: Vec3, a: Vec3) -> (f64, f64) {
parts.iter().flat_map(|(_, trim)| trim.chords.iter()).flat_map(|(x, y)| [*x, *y]).fold((f64::MAX, f64::MIN), |(lo, hi), q| {
let t = (q - p).dot(a);
(lo.min(t), hi.max(t))
})
}
/// Silhouettes of a set of faces on one surface (`parts`: each face's chart and trim, all
/// recognised as `an`): traced once and kept where they lie on any of the faces, so a
/// silhouette along a seam between them is one curve.
pub(crate) fn silhouettes(parts: &[(&S3, &Trim)], an: Analytic, d: Vec3, out: &mut Vec<Sil>) {
let Some(&(s, trim)) = parts.first() else { return };
let mut push = |curve: C3, acc: CurveAccuracy| {
// A traced curve's chords stand up to their deviation off the surface.
let slack = match acc {
CurveAccuracy::Sampled { max_sample_deviation_mm: dev } => dev,
_ => 0.,
};
let mut ranges: Vec<[f64; 2]> = parts.iter().flat_map(|(s, trim)| inside_ranges(&curve, s, trim, slack)).collect();
ranges.sort_by(|a, b| a[0].total_cmp(&b[0]));
let [t0, t1] = curve.domain();
let mut merged: Vec<[f64; 2]> = Vec::with_capacity(ranges.len());
for r in ranges {
match merged.last_mut() {
Some(last) if r[0] <= last[1] + 1e-12 * (t1 - t0) => last[1] = last[1].max(r[1]),
_ => merged.push(r),
}
}
if !merged.is_empty() {
out.push(Sil { curve, acc, ranges: merged });
}
};
match an {
Analytic::Plane => {}
Analytic::Cylinder { p, a, r } => {
let Some(m) = a.cross(d).normalized() else { return };
let (t0, t1) = extent(parts, p, a);
if !(t1 > t0) {
return;
}
for side in [1., -1.] {
let base = p + m * (r * side);
push(C3::Line { a: base + a * t0, b: base + a * t1 }, CurveAccuracy::Analytic);
}
}
Analytic::Cone { apex, a, half } => {
let e1 = perp(a);
let e2 = a.cross(e1);
let (da, d1, d2) = (d.dot(a), d.dot(e1), d.dot(e2));
let rho = (d1 * d1 + d2 * d2).sqrt();
if rho <= 1e-12 {
return;
}
let c = half.tan() * da / rho;
if c.abs() > 1. {
return;
}
let phid = d2.atan2(d1);
let (h0, h1) = extent(parts, apex, a);
let (t0, t1) = (h0.max(0.) / half.cos(), h1.max(0.) / half.cos());
if !(t1 > t0) {
return;
}
for sign in [1., -1.] {
let phi = phid + sign * c.acos();
let g = a * half.cos() + (e1 * phi.cos() + e2 * phi.sin()) * half.sin();
push(C3::Line { a: apex + g * t0, b: apex + g * t1 }, CurveAccuracy::Analytic);
}
}
Analytic::Sphere { c, r } => {
let x = perp(d);
let y = d.cross(x);
push(C3::Arc { o: c, x, y, r: [r, r], range: [0., TAU] }, CurveAccuracy::Analytic);
}
Analytic::Torus { .. } | Analytic::Other => {
// Several faces of one torus may each hold part of its silhouette: trace it on the
// whole torus then. One face keeps its own chart (finer along a narrow strip).
let whole = if parts.len() > 1 { S3::from_analytic(&an, [0., 0.], Vec3::ZERO) } else { None };
let s = whole.as_ref().unwrap_or(s);
let traced = if matches!(s, S3::Rev(_)) { s.rev_silhouette(d) } else { contour(s, trim, d) };
for (curve, dev) in traced {
let acc = if dev == 0. && matches!(curve, C3::Arc { .. }) {
CurveAccuracy::Analytic
} else {
CurveAccuracy::Sampled { max_sample_deviation_mm: dev }
};
push(curve, acc);
}
}
}
}
/// How far a face's surface may lie from the surface its silhouette set is traced on: rounding,
/// 1e-13 of the face's boundary size plus its coordinates (faces join a set when their
/// recognised surfaces agree to 1e-13). Only that close does a ray along the set's silhouette
/// leave the face within the own-face margin (1e-6 of the size): a surface `e` away is crossed
/// about `sqrt(2 R e)` along the ray.
pub(crate) fn join_tol(trim: &Trim) -> f64 {
1e-13 * (1. + (trim.hi - trim.lo).len()) + 1e-15 * trim.lo.len().max(trim.hi.len())
}
/// Parameter ranges of `c` whose points are inside the face, transitions found by bisection.
/// `slack` is how far the curve itself may stand off the surface (a traced curve's chords).
pub(crate) fn inside_ranges(c: &C3, s: &S3, trim: &Trim, slack: f64) -> Vec<[f64; 2]> {
let [t0, t1] = c.domain();
let n = 256;
let mut hint = None;
// The face is closed: a silhouette along its boundary (a seam between two faces of one
// surface, which is not drawn as an edge) is inside. Off the patch is outside: the inversion
// clamps to the chart's domain.
let tol = 2. * trim.chord_error + 1e-9 * (1. + (trim.hi - trim.lo).len());
let on_patch = tol + 2. * join_tol(trim) + 2. * slack;
let state = |t: f64, hint: &mut Option<[f64; 2]>| -> bool {
let p = c.at(t);
let uv = s.inverse(p, *hint);
*hint = uv;
let Some(uv) = uv else { return false };
match trim.contains(uv) {
_ if s.point(uv).dist(p) > on_patch => false,
Some(true) => true,
Some(false) => trim.near(p, tol),
None => false,
}
};
// Uniform samples, plus the parameters nearest to the face's boundary points: an in-face
// interval is bounded by boundary crossings, so a narrow one between two uniform samples
// still gets samples on both of its sides.
let mut ts: Vec<f64> = (0..=n).map(|i| t0 + (t1 - t0) * i as f64 / n as f64).collect();
let near = 0.02 * (trim.hi - trim.lo).len();
ts.extend(trim.chords.iter().step_by(4).filter_map(|(q, _)| closest_param(c, *q).filter(|t| c.at(*t).dist(*q) <= near)));
ts.retain(|t| *t >= t0 && *t <= t1);
ts.sort_by(f64::total_cmp);
ts.dedup_by(|a, b| (*a - *b).abs() <= 1e-12 * (t1 - t0));
let mut hints = Vec::with_capacity(ts.len());
let states: Vec<bool> = ts
.iter()
.map(|t| {
let st = state(*t, &mut hint);
hints.push(hint);
st
})
.collect();
let mut out = Vec::new();
let mut start: Option<f64> = None;
for i in 0..ts.len() {
let (Some(&t), Some(&inside)) = (ts.get(i), states.get(i)) else { continue };
let prev = i.checked_sub(1).and_then(|j| states.get(j).copied());
if prev.is_some_and(|p| p != inside) {
let (mut lo, mut hi) = (ts.get(i - 1).copied().unwrap_or(t), t);
let near = hints.get(i).copied().flatten();
for _ in 0..50 {
let mid = 0.5 * (lo + hi);
let mut h = near;
if state(mid, &mut h) == inside {
hi = mid;
} else {
lo = mid;
}
}
let edge = 0.5 * (lo + hi);
if inside {
start = Some(edge);
} else if let Some(s0) = start.take() {
out.push([s0, edge]);
}
} else if i == 0 && inside {
start = Some(t);
}
}
if let Some(s0) = start {
out.push([s0, t1]);
}
out.retain(|r| r[1] - r[0] > 1e-12 * (t1 - t0));
out
}
/// The parameter of the point of `c` nearest to `p` (lines and arcs exactly, polylines at a
/// vertex), when `p` projects inside the curve's domain.
fn closest_param(c: &C3, p: Vec3) -> Option<f64> {
match c {
C3::Line { a, b } => {
let e = *b - *a;
let l2 = e.len2();
(l2 > 0.).then(|| ((p - *a).dot(e) / l2).clamp(0., 1.))
}
C3::Arc { o, x, y, r, range } => {
let w = p - *o;
let t = (w.dot(*y) / r[1]).atan2(w.dot(*x) / r[0]);
// The turn of t that falls in the range, if any.
let k = ((range[0] - t) / TAU).ceil();
let t = t + k * TAU;
(t <= range[1]).then_some(t)
}
// The nearest point on any chord, not just the nearest vertex: a face narrower than
// the chords still gets parameters inside it.
C3::Polyline(pts) => pts
.windows(2)
.enumerate()
.filter_map(|(i, w)| {
let (a, b) = (*w.first()?, *w.get(1)?);
let e = b - a;
let f = if e.len2() > 0. { ((p - a).dot(e) / e.len2()).clamp(0., 1.) } else { 0. };
Some((i as f64 + f, (a + e * f).dist(p)))
})
.min_by(|x, y| x.1.total_cmp(&y.1))
.map(|(t, _)| t),
C3::Spline { .. } => None,
}
}
/// Zero set of `normal . d` on a chart grid over the face's chart box.
fn contour(s: &S3, trim: &Trim, d: Vec3) -> Vec<(C3, f64)> {
let Some((lo, hi)) = trim.uv_bounds() else { return Vec::new() };
let dom = s.domain();
let lo = [lo[0].max(dom[0][0]), lo[1].max(dom[1][0])];
let hi = [hi[0].min(dom[0][1]), hi[1].min(dom[1][1])];
if !(hi[0] > lo[0] && hi[1] > lo[1]) {
return Vec::new();
}
let g = 128usize;
let uv = |i: usize, j: usize| [lo[0] + (hi[0] - lo[0]) * i as f64 / g as f64, lo[1] + (hi[1] - lo[1]) * j as f64 / g as f64];
let f = |p: [f64; 2]| s.normal(p).map_or(f64::NAN, |n| n.dot(d));
let vals: Vec<f64> = (0..=g).flat_map(|j| (0..=g).map(move |i| (i, j))).map(|(i, j)| f(uv(i, j))).collect();
let val = |i: usize, j: usize| vals.get(j * (g + 1) + i).copied().unwrap_or(f64::NAN);
let zero = |a: [f64; 2], b: [f64; 2], fa: f64| -> [f64; 2] {
let (mut lo, mut hi, mut flo) = (a, b, fa);
for _ in 0..50 {
let m = [0.5 * (lo[0] + hi[0]), 0.5 * (lo[1] + hi[1])];
let fm = f(m);
if fm.is_nan() {
break;
}
if (fm > 0.) == (flo > 0.) {
lo = m;
flo = fm;
} else {
hi = m;
}
}
[0.5 * (lo[0] + hi[0]), 0.5 * (lo[1] + hi[1])]
};
// Edge ids: horizontal (i,j)-(i+1,j) = 2k, vertical (i,j)-(i,j+1) = 2k+1, k = j*(g+1)+i.
let mut points = std::collections::BTreeMap::new();
let mut crossing = |id: usize, a: (usize, usize), b: (usize, usize)| -> Option<usize> {
let (fa, fb) = (val(a.0, a.1), val(b.0, b.1));
if fa.is_nan() || fb.is_nan() || (fa > 0.) == (fb > 0.) {
return None;
}
points.entry(id).or_insert_with(|| zero(uv(a.0, a.1), uv(b.0, b.1), fa));
Some(id)
};
let mut links: std::collections::BTreeMap<usize, Vec<usize>> = std::collections::BTreeMap::new();
for j in 0..g {
for i in 0..g {
let k = |i: usize, j: usize| j * (g + 1) + i;
let e = [
crossing(2 * k(i, j), (i, j), (i + 1, j)),
crossing(2 * k(i + 1, j) + 1, (i + 1, j), (i + 1, j + 1)),
crossing(2 * k(i, j + 1), (i, j + 1), (i + 1, j + 1)),
crossing(2 * k(i, j) + 1, (i, j), (i, j + 1)),
];
let found: Vec<usize> = e.iter().flatten().copied().collect();
let pairs: Vec<(usize, usize)> = match found.as_slice() {
[a, b] => vec![(*a, *b)],
[a, b, c, dd] => {
let centre = f([0.5 * (uv(i, j)[0] + uv(i + 1, j)[0]), 0.5 * (uv(i, j)[1] + uv(i, j + 1)[1])]);
if (centre > 0.) == (val(i, j) > 0.) { vec![(*a, *b), (*c, *dd)] } else { vec![(*a, *dd), (*b, *c)] }
}
_ => Vec::new(),
};
for (a, b) in pairs {
links.entry(a).or_default().push(b);
links.entry(b).or_default().push(a);
}
}
}
// Chain the links into polylines.
let mut used = std::collections::BTreeSet::new();
let mut chains = Vec::new();
let starts: Vec<usize> = links.iter().filter(|(_, v)| v.len() == 1).map(|(k, _)| *k).chain(links.keys().copied()).collect();
for s0 in starts {
if used.contains(&s0) {
continue;
}
let mut chain = vec![s0];
used.insert(s0);
let mut cur = s0;
while let Some(next) = links.get(&cur).and_then(|v| v.iter().find(|n| !used.contains(*n)).copied()) {
used.insert(next);
chain.push(next);
cur = next;
}
// Close loops back to their start.
if chain.len() > 2 && links.get(&cur).is_some_and(|v| v.contains(&s0)) {
chain.push(s0);
}
if chain.len() >= 2 {
chains.push(chain);
}
}
let mut out = Vec::new();
for chain in chains {
let uvs: Vec<[f64; 2]> = chain.iter().filter_map(|id| points.get(id).copied()).collect();
let pts: Vec<Vec3> = uvs.iter().map(|p| s.point(*p)).collect();
// Deviation: refine each chord's midpoint onto the curve across the chord.
let mut dev = 0f64;
for (w, q) in uvs.windows(2).zip(pts.windows(2)) {
let (Some(a), Some(b), Some(pa), Some(pb)) = (w.first(), w.get(1), q.first(), q.get(1)) else { continue };
let m = [0.5 * (a[0] + b[0]), 0.5 * (a[1] + b[1])];
let n = [-(b[1] - a[1]), b[0] - a[0]];
let (x, y) = ([m[0] - n[0], m[1] - n[1]], [m[0] + n[0], m[1] + n[1]]);
let (fx, fy) = (f(x), f(y));
if fx.is_nan() || fy.is_nan() || (fx > 0.) == (fy > 0.) {
continue;
}
let exact = s.point(zero(x, y, fx));
dev = dev.max(exact.dist_to_segment(*pa, *pb));
}
out.push((C3::Polyline(pts), dev));
}
out
}
crates/geom/src/hlr/surface.rs
+1298
View File
File diff suppressed because it is too large Load Diff
crates/geom/src/hlr/trim.rs
+382
View File
@@ -0,0 +1,382 @@
//! Face trimming in the surface chart. Boundary loops are sampled from the exact edge curves,
//! unwrapped across periodic directions, and a point is classified by the parity of boundary
//! crossings on a path to a point known to be inside. Seams (an edge used twice by one face)
//! are crossed twice, so they cancel without special cases.
use super::curve::C3;
use super::surface::S3;
use crate::Vec3;
use crate::brep::Face;
#[derive(Clone, Debug)]
pub(crate) struct Trim {
/// Unwrapped chart segments of all loops.
segs: Vec<([f64; 2], [f64; 2])>,
period: [Option<f64>; 2],
/// A point inside the face; `None` when the face covers its whole chart.
inside: Option<[f64; 2]>,
/// 3D boundary chords, for distance-to-boundary margins.
pub chords: Vec<(Vec3, Vec3)>,
/// Largest distance between an exact boundary curve and its chords.
pub chord_error: f64,
pub lo: Vec3,
pub hi: Vec3,
/// Why classification is unavailable, if it is.
pub unsupported: Option<&'static str>,
/// Index for the two-leg parity path q -> (r.u, q.v) -> r.
index: Option<Index>,
/// `chords` on a grid, for boundaries with many of them.
grid: Option<ChordGrid>,
}
/// Boundary chords bucketed on a uniform grid over the boundary's box, so a point is checked
/// against the chords around it, not all of them (a circle alone has `trim_samples` per break).
/// Cells are at least as long as the longest chord, so a chord lies in at most two per axis.
#[derive(Clone, Debug)]
struct ChordGrid {
lo: Vec3,
cell: f64,
n: [usize; 3],
/// (cell, chord) pairs sorted by cell.
entries: Vec<(u32, u32)>,
}
impl ChordGrid {
fn new(chords: &[(Vec3, Vec3)], lo: Vec3, hi: Vec3) -> Option<ChordGrid> {
// A few chords are checked faster one by one than through a grid.
if chords.len() < 64 || u32::try_from(chords.len()).is_err() {
return None;
}
let ext = hi - lo;
let longest = chords.iter().map(|(a, b)| a.dist(*b)).fold(0., f64::max);
let cell = longest.max(ext.len() / 256.).max(ext.x.max(ext.y).max(ext.z) / 512.);
if !(cell > 0. && cell.is_finite()) {
return None;
}
let n = [ext.x, ext.y, ext.z].map(|e| ((e / cell).floor() as usize).saturating_add(1).min(513));
let mut grid = ChordGrid { lo, cell, n, entries: Vec::with_capacity(chords.len() * 2) };
for (i, (a, b)) in (0u32..).zip(chords) {
let (c0, c1) = (grid.cell_of(a.min(*b)), grid.cell_of(a.max(*b)));
for x in c0[0]..=c1[0] {
for y in c0[1]..=c1[1] {
for z in c0[2]..=c1[2] {
grid.entries.push((grid.key([x, y, z]), i));
}
}
}
}
grid.entries.sort_unstable();
Some(grid)
}
fn cell_of(&self, p: Vec3) -> [usize; 3] {
let at = |v: f64, l: f64, n: usize| (((v - l) / self.cell).floor().max(0.) as usize).min(n.saturating_sub(1));
[at(p.x, self.lo.x, self.n[0]), at(p.y, self.lo.y, self.n[1]), at(p.z, self.lo.z, self.n[2])]
}
fn key(&self, c: [usize; 3]) -> u32 {
// At most 513 cells per axis: the key fits in 28 bits.
((c[0] * self.n[1] + c[1]) * self.n[2] + c[2]) as u32
}
/// Whether a chord touching a cell of the box `p` +- `tol` is within `tol` of `p`; `None`
/// when that box covers too many cells to be worth it.
fn near(&self, chords: &[(Vec3, Vec3)], p: Vec3, tol: f64) -> Option<bool> {
// A little wider than `tol`, so rounding in the distance cannot pick a chord whose
// nearest point lies in a cell left out.
let w = tol * (1. + 1e-6);
let t = Vec3::new(w, w, w);
let (c0, c1) = (self.cell_of(p - t), self.cell_of(p + t));
let span = |k: usize| c1[k].saturating_sub(c0[k]).saturating_add(1);
if span(0).saturating_mul(span(1)).saturating_mul(span(2)) > 64 {
return None;
}
for x in c0[0]..=c1[0] {
for y in c0[1]..=c1[1] {
for z in c0[2]..=c1[2] {
let k = self.key([x, y, z]);
let first = self.entries.partition_point(|e| e.0 < k);
let close = self
.entries
.get(first..)
.unwrap_or_default()
.iter()
.take_while(|e| e.0 == k)
.any(|(_, i)| chords.get(*i as usize).is_some_and(|(a, b)| p.dist_to_segment(*a, *b) <= tol));
if close {
return Some(true);
}
}
}
}
Some(false)
}
}
/// Segment copies shifted across periods near the inside point, the sorted crossings of the
/// vertical line u = r.u, and v-buckets of the copies for the horizontal leg.
#[derive(Clone, Debug)]
struct Index {
copies: Vec<([f64; 2], [f64; 2])>,
vcross: Vec<f64>,
buckets: Vec<Vec<u32>>,
v0: f64,
dv: f64,
}
fn spans(x: f64, a: f64, b: f64) -> bool {
(a <= x && x < b) || (b <= x && x < a)
}
fn index(segs: &[([f64; 2], [f64; 2])], period: [Option<f64>; 2], r: [f64; 2]) -> Index {
let shifts = |k: usize| -> Vec<f64> {
let Some(p) = period[k] else { return vec![0.] };
let (lo, hi) = segs.iter().fold((f64::MAX, f64::MIN), |(lo, hi), (a, b)| (lo.min(a[k]).min(b[k]), hi.max(a[k]).max(b[k])));
let k0 = ((r[k] - p - hi) / p).floor() as i64;
let k1 = ((r[k] + p - lo) / p).ceil() as i64;
(k0.max(-8)..=k1.min(8)).map(|i| i as f64 * p).collect()
};
let mut copies = Vec::new();
for dx in shifts(0) {
for dy in shifts(1) {
copies.extend(segs.iter().map(|(a, b)| ([a[0] + dx, a[1] + dy], [b[0] + dx, b[1] + dy])));
}
}
let mut vcross: Vec<f64> =
copies.iter().filter(|(a, b)| spans(r[0], a[0], b[0])).map(|(a, b)| a[1] + (b[1] - a[1]) * (r[0] - a[0]) / (b[0] - a[0])).collect();
vcross.sort_by(f64::total_cmp);
let (v0, v1) = copies.iter().fold((f64::MAX, f64::MIN), |(lo, hi), (a, b)| (lo.min(a[1]).min(b[1]), hi.max(a[1]).max(b[1])));
let n = ((copies.len() as f64).sqrt().ceil() as usize).clamp(1, 4096);
let dv = ((v1 - v0) / n as f64).max(1e-300);
let mut buckets = vec![Vec::new(); n];
for (i, (a, b)) in copies.iter().enumerate() {
let k0 = (((a[1].min(b[1]) - v0) / dv).floor().max(0.) as usize).min(n - 1);
let k1 = (((a[1].max(b[1]) - v0) / dv).floor().max(0.) as usize).min(n - 1);
for bucket in buckets.iter_mut().take(k1 + 1).skip(k0) {
bucket.push(i as u32);
}
}
Index { copies, vcross, buckets, v0, dv }
}
fn lift(x: f64, near: f64, p: Option<f64>) -> f64 {
match p {
Some(p) if p > 0. => x + p * ((near - x) / p).round(),
_ => x,
}
}
pub(crate) fn build(face: &Face, edges: &[C3], surface: &S3, per_break: usize) -> Trim {
let period = surface.periods();
let mut segs = Vec::new();
let mut seam = Vec::new();
let mut chords = Vec::new();
let mut chord_error = 0f64;
let (mut lo, mut hi) = (Vec3::new(f64::MAX, f64::MAX, f64::MAX), Vec3::new(f64::MIN, f64::MIN, f64::MIN));
let mut uses = std::collections::BTreeMap::new();
for u in face.loops.iter().flatten() {
*uses.entry(u.edge).or_insert(0usize) += 1;
}
let mut area_u = 0.; // -loop integral of v du
let mut area_v = 0.; // loop integral of u dv
let mut wraps = [false; 2];
let mut unsupported = None;
for lp in &face.loops {
let mut pts: Vec<([f64; 2], bool)> = Vec::new();
let mut hint: Option<[f64; 2]> = None;
for u in lp {
let Some(c) = edges.get(u.edge as usize) else {
unsupported = Some("edge reference");
continue;
};
let is_seam = uses.get(&u.edge).copied().unwrap_or(0) > 1;
let mut ts = Vec::new();
let b = c.breaks();
for w in b.windows(2) {
let (Some(t0), Some(t1)) = (w.first(), w.get(1)) else { continue };
for i in 0..per_break {
ts.push(t0 + (t1 - t0) * i as f64 / per_break as f64);
}
}
if let Some(last) = b.last() {
ts.push(*last);
}
if !u.forward {
ts.reverse();
}
let mut prev3: Option<(f64, Vec3)> = None;
for t in ts {
let p = c.at(t);
lo = lo.min(p);
hi = hi.max(p);
if let Some((tp, q)) = prev3 {
chords.push((q, p));
let m = c.at(0.5 * (t + tp));
chord_error = chord_error.max(m.dist_to_segment(q, p));
}
prev3 = Some((t, p));
let Some(mut uv) = surface.inverse(p, hint) else {
unsupported = Some("boundary point off its surface");
continue;
};
if let Some(h) = hint {
uv = [lift(uv[0], h[0], period[0]), lift(uv[1], h[1], period[1])];
}
hint = Some(uv);
pts.push((uv, is_seam));
}
}
let (Some(first), Some(last)) = (pts.first().copied(), pts.last().copied()) else { continue };
// Close the loop at the lifted copy of its start (a wrapping loop ends a period away).
let close = [lift(first.0[0], last.0[0], period[0]), lift(first.0[1], last.0[1], period[1])];
for k in 0..2 {
if (close[k] - first.0[k]).abs() > 1e-9 {
wraps[k] = true;
}
}
pts.push((close, last.1));
for w in pts.windows(2) {
let (Some((a, _)), Some((b, s))) = (w.first(), w.get(1)) else { continue };
area_u -= 0.5 * (a[1] + b[1]) * (b[0] - a[0]);
area_v += 0.5 * (a[0] + b[0]) * (b[1] - a[1]);
segs.push((*a, *b));
seam.push(*s);
}
}
let area = match wraps {
[true, true] => {
unsupported = unsupported.or(Some("loops wrap in both chart directions"));
0.
}
[false, true] => area_v,
_ => area_u,
};
// Loops that wind once around a periodic direction without coming back (a hemisphere's
// equator) close against the chart's own edge, so their area depends on where the chart
// starts: there the face is taken on the loops' left (the B-rep convention), and below the
// side within the surface's domain wins either way.
let net = segs.iter().fold([0f64; 2], |n, (a, b)| [n[0] + b[0] - a[0], n[1] + b[1] - a[1]]);
let open = (0..2).any(|k| period[k].is_some_and(|p| net[k].abs() > 0.5 * p));
let sigma = if open { 1. } else { area.signum() };
let dom = surface.domain();
// Inside point: beside the longest non-seam segment, on the face's side.
let mut inside = None;
let real: Vec<usize> = (0..segs.len()).filter(|i| !seam.get(*i).copied().unwrap_or(true)).collect();
if !real.is_empty() && (open || area != 0.) {
let len = |i: usize| segs.get(i).map_or(0., |(a, b)| ((b[0] - a[0]).powi(2) + (b[1] - a[1]).powi(2)).sqrt());
let best = real.iter().copied().max_by(|a, b| len(*a).total_cmp(&len(*b)));
if let Some((a, b)) = best.and_then(|i| segs.get(i)) {
let m = [0.5 * (a[0] + b[0]), 0.5 * (a[1] + b[1])];
let l = ((b[0] - a[0]).powi(2) + (b[1] - a[1]).powi(2)).sqrt();
let n = [-(b[1] - a[1]) / l * sigma, (b[0] - a[0]) / l * sigma];
// Stay well inside: a quarter of the distance to any other boundary piece.
let mut room = l;
for (j, (c, e)) in segs.iter().enumerate() {
if Some(j) == best {
continue;
}
for (dx, dy) in shifts(period) {
let (c2, e2) = ([c[0] + dx, c[1] + dy], [e[0] + dx, e[1] + dy]);
room = room.min(point_segment(m, c2, e2));
}
}
let delta = 0.25 * room;
// At most half way to the chart's domain edge, so a face hugging that edge (a cone's
// tip) keeps its point; a loop on the edge itself has no room on the outer side.
let step = |dir: [f64; 2]| {
let mut s = delta;
for (((d, &c), &mk), p) in dom.iter().zip(&dir).zip(&m).zip(&period) {
if p.is_some() {
continue;
}
if c < 0. {
s = s.min(0.5 * (mk - d[0]) / -c);
} else if c > 0. {
s = s.min(0.5 * (d[1] - mk) / c);
}
}
s
};
let flip = [-n[0], -n[1]];
let pick = [(n, step(n)), (flip, step(flip))].into_iter().find(|(_, s)| *s > 0. && s.is_finite());
if let Some((dir, s)) = pick {
inside = Some([m[0] + dir[0] * s, m[1] + dir[1] * s]);
}
}
if inside.is_none() {
unsupported = unsupported.or(Some("no room for an inside reference point"));
}
}
let index = inside.map(|r| index(&segs, period, r));
let grid = ChordGrid::new(&chords, lo, hi);
Trim { segs, period, inside, chords, chord_error, lo, hi, unsupported, index, grid }
}
fn shifts(period: [Option<f64>; 2]) -> Vec<(f64, f64)> {
let ks = |p: Option<f64>| p.map_or(vec![0.], |p| vec![-p, 0., p]);
let mut out = Vec::new();
for dx in ks(period[0]) {
for dy in ks(period[1]) {
out.push((dx, dy));
}
}
out
}
fn point_segment(p: [f64; 2], a: [f64; 2], b: [f64; 2]) -> f64 {
let (ex, ey) = (b[0] - a[0], b[1] - a[1]);
let l2 = ex * ex + ey * ey;
let t = if l2 > 0. { (((p[0] - a[0]) * ex + (p[1] - a[1]) * ey) / l2).clamp(0., 1.) } else { 0. };
((a[0] + t * ex - p[0]).powi(2) + (a[1] + t * ey - p[1]).powi(2)).sqrt()
}
impl Trim {
/// Whether chart point `uv` lies inside the face; `None` when the trim is unsupported.
/// Parity of boundary crossings on the path q -> (r.u, q.v) -> r to the inside point r.
pub fn contains(&self, uv: [f64; 2]) -> Option<bool> {
if self.unsupported.is_some() {
return None;
}
let (Some(r), Some(ix)) = (self.inside, self.index.as_ref()) else { return Some(true) };
let q = [lift(uv[0], r[0], self.period[0]), lift(uv[1], r[1], self.period[1])];
let mut crossings = 0usize;
// Horizontal leg at v = q.v between q.u and r.u.
let (ulo, uhi) = (q[0].min(r[0]), q[0].max(r[0]));
let k = (((q[1] - ix.v0) / ix.dv).floor().max(0.) as usize).min(ix.buckets.len().saturating_sub(1));
for i in ix.buckets.get(k).map(Vec::as_slice).unwrap_or(&[]) {
let Some((a, b)) = ix.copies.get(*i as usize) else { continue };
if spans(q[1], a[1], b[1]) {
let u = a[0] + (b[0] - a[0]) * (q[1] - a[1]) / (b[1] - a[1]);
if u > ulo && u <= uhi {
crossings += 1;
}
}
}
// Vertical leg at u = r.u between q.v and r.v.
let (vlo, vhi) = (q[1].min(r[1]), q[1].max(r[1]));
let first = ix.vcross.partition_point(|v| *v <= vlo);
let last = ix.vcross.partition_point(|v| *v <= vhi);
crossings += last.saturating_sub(first);
Some(crossings.is_multiple_of(2))
}
/// The chart box of the unwrapped boundary.
pub fn uv_bounds(&self) -> Option<([f64; 2], [f64; 2])> {
let mut it = self.segs.iter().flat_map(|(a, b)| [a, b]);
let first = *it.next()?;
Some(it.fold((first, first), |(lo, hi), p| ([lo[0].min(p[0]), lo[1].min(p[1])], [hi[0].max(p[0]), hi[1].max(p[1])])))
}
/// Whether a 3D point is within `tol` of the sampled boundary: the boundary's box first,
/// then the chords around the point (all of them for a short boundary or a wide `tol`).
pub fn near(&self, p: Vec3, tol: f64) -> bool {
let out = (self.lo - p).max(p - self.hi).max(Vec3::new(0., 0., 0.));
if !(out.len() <= tol) {
return false;
}
match self.grid.as_ref().and_then(|g| g.near(&self.chords, p, tol)) {
Some(close) => close,
None => self.chords.iter().any(|(a, b)| p.dist_to_segment(*a, *b) <= tol),
}
}
}
crates/geom/src/hlr_tests.rs
+2213
View File
File diff suppressed because it is too large Load Diff
crates/geom/src/lib.rs
+1
View File
@@ -24,6 +24,7 @@ pub mod env;
pub mod fill;
pub mod fill_general;
mod grid;
pub mod hlr;
mod mesh;
pub mod mesh_distance;
pub mod network;
docs/drawing-hlr.md
+198
View File
@@ -0,0 +1,198 @@
# Hidden lines (W06a)
This is the geometry core of Drawing package W06 from the [published plan][plan] in
[#860][rfc]: [#2455][issue] / [PR #2457][pr]. It turns the subset that the W05a
investigation ([#2105][w05-issue], [draft #2106][w05-pr], [findings][findings])
showed to be ready into code: exact orthographic hidden lines for placed bodies whose
faces are planes, cylinders, cones, spheres or tori.
The code is `crates/geom/src/hlr` (layer L0). It reads the neutral `geom::brep`
records, uses no kernel or UI types, returns errors instead of panicking and builds for
`wasm32-unknown-unknown`. It does not register commands, create views, extract bodies
from the kernel or draw sections. The read-only kernel adapter follows the W02
extraction contract; until then callers build `Brep` records themselves.
## Using it
```text
prepare(&Brep, &Options) -> Result<Prepared, HlrError> once per body
hidden_lines(&[Prepared], &[Placement], &View, &Options, &cancel) -> Result<Drawing, HlrError>
```
- `prepare` checks the body (`Brep::validate`), recognises each face's surface, builds
its trim and classifies each edge as sharp, tangent or seam. `Prepared::report` says
what it found: surfaces and curves by kind, tangent and seam edges, faces whose trim is
unsupported, freeform faces and the largest boundary chord error.
- A `Placement` puts a prepared body in the scene: `world = rotation * local +
translation`, the rotation given by its columns. Several placements may share one
prepared body. A placement that is not a proper rotation (a mirror, a scale) is refused
with `HlrError::Invalid`.
- `View::new(direction_of_sight, up)` is an orthographic view. Screen coordinates are
`(p . x, p . y)` in millimetres, with `x = direction x up` normalised and `y = x x
direction`.
- `cancel` is polled during the work; when it returns true the call ends with
`HlrError::Cancelled`. `prepare_cancellable(&Brep, &Options, &cancel)` is `prepare` with
such a callback, polled between edges, faces and surface comparisons.
The `Drawing` holds one `Line` per drawable edge and per silhouette curve of each
placement, in order of placement. A line names its `placement` and `source` (`Edge { edge, tangent }`
or `Silhouette { face }`) and is split into `pieces`. Each `Piece` has:
| Field | Meaning |
| --- | --- |
| `t` | Parameter interval on the source curve. |
| `visibility` | `Visible`, `Hidden` or `Uncertain`. |
| `covered` | Another piece already draws the same projected curve here; draw only uncovered pieces. |
| `geometry` | The piece in screen coordinates as `Seg2` segments. |
| `accuracy` | `CurveAccuracy` of `geometry` (the [W01a vocabulary](drawing-geometry.md)). |
`Drawing::unsupported` lists the faces whose silhouettes were not drawn and why;
`Drawing::stats` counts curves, crossings, pieces, rays and uncertain pieces, with phase
timings where the platform has a clock (zero on `wasm32-unknown-unknown`).
## What is exact
| Source | Geometry | Accuracy |
| --- | --- | --- |
| Line edges; cylinder and cone silhouette rulings | one `Seg2::Line` | `Analytic` |
| Circle and ellipse edges; sphere silhouettes; torus silhouettes seen along the axis | `Seg2::Conic` arcs of at most 90 degrees | `Analytic` |
| Spline and intersection edges | chords between exact points, halved until each piece's control hull is within the sag target | `Sampled` (a bound from that hull) |
| Torus silhouettes in other views | chords between traced points, halved to 1e-7 of the torus's size | `Sampled` (measured deviation) |
An orthographic projection is affine, so a circle projects to an ellipse, which a
rational quadratic arc represents exactly. Chords are refined until their sag on screen
is at most 1e-7 of the scene size. A fitted edge record also carries its own bound, and
the visibility work allows for it.
Surfaces are recognised whether they are stored as analytic records, as revolutions of
lines and circles or as spline patches (to 1e-9 relative; a patch whose control points lie
in one plane is that plane). A spline curve or patch is taken for a circle, quadric or torus
only when it lies on it at enough points of every knot span to make the check exact for its
degree, with its largest distance there, scaled by how far apart its weights are, within
the tolerance: weights near zero cannot hide a swerve, and widely spread weights mean chords.
A recognised patch is trimmed and intersected in its own analytic chart. Faces that lie on one surface (a sphere
stored as two hemispheres, a cylinder split in two) share their silhouettes: the edges
between them are seams and are not drawn, and a silhouette running along such a seam is
one curve. A face joins such a set only when the first face's surface passes through its
own points to rounding (1e-13 of its size); recognition alone groups surfaces more loosely,
and such faces keep silhouettes of their own.
## Visibility
Projected curves are split where they cross or touch, found on a screen grid and
refined on the exact curves. Each piece is classified by rays toward the viewer against
every trimmed face: closed forms for planes and quadrics, a sampled quartic for tori.
A ray that grazes a face, meets one within its boundary tolerance or meets another body
right at its start makes the piece `Uncertain`; nothing is guessed.
Coincident curves (an edge and a silhouette on top of each other, two touching bodies)
follow the W05a policy: visible wins and duplicates are marked `covered`. The tolerance
comes from the curves' own accuracy (chord sag and fit bounds) plus a rounding allowance
of 1e-7 mm and 1e-9 of the scene size, never from the scene size alone, so small parts in
large scenes keep their detail.
`prepare` moves each body to its own centre, and the work is done in a frame centred on
the scene and shifted back at the end, so coordinates stay small: a part 1e7 mm from the
origin, whether placed there or modelled there, is classified as at the origin.
## Freeform faces
Faces that are not planes, cylinders, cones, spheres or tori are freeform. By default
they are listed in `Drawing::unsupported` with the reason `freeform face`, their
silhouettes are not drawn, and a piece whose ray meets one is `Uncertain`, never
`Hidden`. Without its silhouettes nothing marks where such a face's outline crosses a curve,
so curves are cut at the face's screen box and every piece inside the box that can reach
behind its nearest point (bounded by the curve's control points, not by samples) is
`Uncertain` too, the face's own edges included (faces left out for an unsupported trim
likewise). With `Options::freeform` set their silhouettes are contoured and their ray
tests use seeded Newton iterations: slower, `Sampled`, and subject to the misses listed
in the W05a findings.
## Limits
- `Options::max_tracks` (curves) and `Options::max_rays` bound the work, as do internal
caps: crossing grid entries (50 million), crossing candidate pairs (500 million, counted
before the search, so no set of pairs met grows with them), crossings (50 million) and
the surface comparisons that group a body's faces (50 million). Exceeding one returns
`HlrError::Budget`.
- `Options::trim_samples` is 1 to `MAX_TRIM_SAMPLES` (4096), and every coordinate of a
placed body at most `MAX_COORDINATE` (1e12 mm, the bound `geom::section` and
`geom::brep` put on stored numbers); otherwise the call returns `HlrError::Invalid`.
- `Options::tangent_edges` (default on) draws edges between tangent faces, flagged
`tangent: true`; drafting policy for them is W17.
- Tolerances are relative to the faces and the scene. They are meant for faces up to
about 10 m across and features down to about 1e-3 mm. A spline patch only a sliver wide
next to its radii (a few millionths of them) cannot be told from its neighbours' surfaces
reliably and is reported as freeform. Distinct surfaces that
recognition takes for one (closer than 1e-9 of their size or position, a nearly flat
cone's far apex included) share the own-face margin of 1e-6 of the face's size when a
ray leaves one of them.
- Not in this slice: parallel per-curve work and the shared work limits of
`geom::section`, the kernel adapter (after W02), perspective, section regions (W03),
views and commands (W07), and mirrored placements. Of the records the kernel now saves,
planes through three points (`PlaneThrough`) are drawn as planes; a body with a procedural
`SurfaceIntersection` curve or a `PlacedRevolution` surface (whose placement may scale it
unevenly) is refused by `prepare` with `HlrError::Unsupported` until the adapter handles
them.
## Fixtures
`crates/geom/src/hlr_tests.rs` builds bodies from records and checks analytic positions
and intervals:
- a cube in an isometric view hides its three back edges; a box in front hides the
middle of an edge over exactly the covered interval;
- a cylinder's silhouette rulings are exact lines, including one along the seam; a cone
frustum seen from the front has exact rulings and draws each edge-on circle once;
- a sphere's outline is exact conic arcs at its radius in three views;
- a disc bounded by a rational circle with unevenly spaced knots, with a knot span narrower
than rounding, or with a parameterisation that races through its first turn is drawn as
arcs turning the full 360 degrees; splines off their circle between evenly spaced samples
or inside such a narrow span are drawn as chords; exact arcs that turn almost all of
their 358 degrees in the first hundredth of their parameter, or several radians in a
knot span of 1e-6 or 1e-11, keep the turn of every piece;
- a torus stored as two halves, seen face on, is drawn at both equators, one silhouette
each; an edge under it is hidden over exactly the intervals the tube covers; seen
obliquely its outline spans R + r on both sides with a deviation under 1e-3 mm;
- cone frusta and narrow cone sectors a hair off one axis or with axes 3e-7 rad apart, and
spherical caps through one circle with centres 5e-7 mm apart, keep their own outlines;
- a cylinder 1e5 and 1e7 mm from the origin draws as at the origin; a cone trimmed to its
tip 1 mm below the apex hides nothing above it;
- a freeform patch is reported and makes the edges it covers uncertain, or hidden when
freeform faces are enabled; a flat spline patch hides them for certain; a long edge behind
a freeform tube is uncertain where the tube may cover it and visible elsewhere, and so is
a sloping edge that only dips behind a freeform patch near its end, while one that leaves
the patch's box below its level but crosses the box above it stays visible;
- a curve that swerves 2.9 mm off its circle where its weight is 1e-12, and a nearly straight
cubic 0.002 mm off a circle of radius 1e7, are drawn as chords; an edge whose depth dips
behind a freeform patch between its samples is uncertain there, and so is an arc seen
along its plane whose projection folds back under a freeform patch between two samples,
and the edge of a freeform sheet that curls back over it 0.005 mm above;
- a torus patch stored as a rational spline is recognised and draws its outer equator as an
exact arc, and one with a bump in a knot span 0.002 of its parameters wide is not taken for
a torus; a big circle passing under a thin freeform tube in a coarse scene is uncertain
there;
- two strips of a torus 0.0003 rad wide keep their 0.006 mm of outline;
- two touching boxes draw each shared edge once;
- a box whose top plane record stands one rounding (or 1e-11 of 100 mm) off its edges still
hides its back edges, and so does a 0.001 mm cube 1e7 mm out whose top plane is one step
above its edges; a plate's sixteen holes each hide the same part of their bottom circle;
eight coaxial cylinders seen along their axis show one circle; 300 boxes apart on screen
are each drawn as the first one; an arc that bulges across a line 1 mm off its chord still
crosses it twice; the trim's chord grid agrees with a scan of every chord; surface groups
agree with comparing every pair of faces, and preparing can be cancelled;
- mirrored placements, bad body indices and bad views are refused; cancellation and
budgets return errors;
- a box stored with planes through three points draws the same bits as one with
axis-stored planes; placed revolutions and procedural intersection curves are refused as
not supported yet;
- torus ray roots for thin, large and nearly tangent cases, and torus silhouettes near
the axis direction.
[rfc]: https://github.com/tangency-labs/solidcraft/issues/860
[plan]: https://github.com/tangency-labs/solidcraft/issues/860#issuecomment-6234779815
[findings]: https://github.com/tangency-labs/solidcraft/issues/860#issuecomment-6243691483
[issue]: https://github.com/tangency-labs/solidcraft/issues/2455
[pr]: https://github.com/tangency-labs/solidcraft/pull/2457
[w05-issue]: https://github.com/tangency-labs/solidcraft/issues/2105
[w05-pr]: https://github.com/tangency-labs/solidcraft/pull/2106