diff --git a/src/modules/draw/draw/hatch.rs b/src/modules/draw/draw/hatch.rs index 05d44e7cc..da673f3fc 100644 --- a/src/modules/draw/draw/hatch.rs +++ b/src/modules/draw/draw/hatch.rs @@ -50,7 +50,7 @@ enum Mode { // ── CPU point-in-polygon (ray casting) ──────────────────────────────────── -fn point_in_polygon(p: [f32; 2], poly: &[[f32; 2]]) -> bool { +fn point_in_polygon(p: [f64; 2], poly: &[[f64; 2]]) -> bool { let n = poly.len(); if n < 3 { return false; @@ -73,25 +73,26 @@ fn point_in_polygon(p: [f32; 2], poly: &[[f32; 2]]) -> bool { /// Shoelace-area magnitude of a polygon. Used to pick the smallest enclosing /// outline when a click falls inside several nested boundaries. -fn polygon_area(poly: &[[f32; 2]]) -> f32 { +fn polygon_area(poly: &[[f64; 2]]) -> f64 { let n = poly.len(); if n < 3 { return 0.0; } - let mut a = 0.0_f64; - let mut j = n - 1; - for i in 0..n { - a += (poly[j][0] as f64) * (poly[i][1] as f64) - - (poly[i][0] as f64) * (poly[j][1] as f64); - j = i; + let origin = poly[0]; + let mut a = 0.0; + for i in 1..n - 1 { + let current = poly[i]; + let next = poly[i + 1]; + a += (current[0] - origin[0]) * (next[1] - origin[1]) + - (next[0] - origin[0]) * (current[1] - origin[1]); } - (a * 0.5).abs() as f32 + (a * 0.5).abs() } /// True when every vertex of `inner` lies inside `outer`. Sufficient to /// recognise a closed hatch outline as nested inside another for the common /// rectangle / closed-polyline case. -fn polygon_contains_polygon(outer: &[[f32; 2]], inner: &[[f32; 2]]) -> bool { +fn polygon_contains_polygon(outer: &[[f64; 2]], inner: &[[f64; 2]]) -> bool { if inner.len() < 3 { return false; } @@ -110,10 +111,10 @@ fn polygon_contains_polygon(outer: &[[f32; 2]], inner: &[[f32; 2]]) -> bool { /// * click inside the innermost shape → hatch just that shape, /// * click in a gap → hatch that ring, with the next level in as holes. fn resolve_hatch_rings( - outlines: &[Vec<[f32; 2]>], - p: [f32; 2], + outlines: &[Vec<[f64; 2]>], + p: [f64; 2], ) -> Option>> { - let mut containing: Vec<(usize, f32)> = outlines + let mut containing: Vec<(usize, f64)> = outlines .iter() .enumerate() .filter(|(_, o)| point_in_polygon(p, o)) @@ -127,8 +128,7 @@ fn resolve_hatch_rings( let outer_idx = containing[0].0; let outer = &outlines[outer_idx]; - let mut rings: Vec> = - vec![outer.iter().map(|&[x, y]| [x as f64, y as f64]).collect()]; + let mut rings = vec![outer.clone()]; for (i, o) in outlines.iter().enumerate() { if i == outer_idx { continue; @@ -149,7 +149,7 @@ fn resolve_hatch_rings( && polygon_contains_polygon(x, o) }); if !has_intermediate { - rings.push(o.iter().map(|&[x, y]| [x as f64, y as f64]).collect()); + rings.push(o.clone()); } } Some(rings) @@ -194,14 +194,14 @@ fn rte_boundary(pts: impl Iterator) -> (Vec<[f32; 2]>, [f64; // ── HATCH command ────────────────────────────────────────────────────────── pub struct HatchCommand { - outlines: Vec>, + outlines: Vec>, mode: Mode, manual_pts: Vec, missed: bool, } impl HatchCommand { - pub fn new(outlines: Vec>) -> Self { + pub fn new(outlines: Vec>) -> Self { Self { outlines, mode: Mode::PickInside, @@ -294,8 +294,7 @@ impl CadCommand for HatchCommand { fn on_point(&mut self, pt: DVec3) -> CmdResult { match &self.mode { Mode::PickInside => { - // Hit-test against the f32 outline catalog (screen-space). - let xy = [pt.x as f32, pt.y as f32]; + let xy = [pt.x, pt.y]; match resolve_hatch_rings(&self.outlines, xy) { Some(rings) => { self.missed = false; @@ -375,7 +374,7 @@ impl CadCommand for HatchCommand { // ── GRADIENT command ─────────────────────────────────────────────────────── pub struct GradientCommand { - outlines: Vec>, + outlines: Vec>, mode: Mode, manual_pts: Vec, missed: bool, @@ -386,7 +385,7 @@ pub struct GradientCommand { } impl GradientCommand { - pub fn new(outlines: Vec>) -> Self { + pub fn new(outlines: Vec>) -> Self { Self { outlines, mode: Mode::PickInside, @@ -486,8 +485,7 @@ impl CadCommand for GradientCommand { fn on_point(&mut self, pt: DVec3) -> CmdResult { match &self.mode { Mode::PickInside => { - // Hit-test against the f32 outline catalog (screen-space). - let xy = [pt.x as f32, pt.y as f32]; + let xy = [pt.x, pt.y]; match resolve_hatch_rings(&self.outlines, xy) { Some(rings) => { self.missed = false; @@ -581,12 +579,12 @@ impl CadCommand for GradientCommand { // ── BOUNDARY command ─────────────────────────────────────────────────────── pub struct BoundaryCommand { - outlines: Vec>, + outlines: Vec>, missed: bool, } impl BoundaryCommand { - pub fn new(outlines: Vec>) -> Self { + pub fn new(outlines: Vec>) -> Self { Self { outlines, missed: false, @@ -609,8 +607,7 @@ impl CadCommand for BoundaryCommand { } fn on_point(&mut self, pt: DVec3) -> CmdResult { - // Hit-test against the f32 outline catalog (screen-space). - let xy = [pt.x as f32, pt.y as f32]; + let xy = [pt.x, pt.y]; match resolve_hatch_rings(&self.outlines, xy) { Some(rings) => { self.missed = false; @@ -656,13 +653,13 @@ inventory::submit!(crate::command::CommandRegistration { names: &["HATCH"] }); mod tests { use super::*; - fn rect(x0: f32, y0: f32, x1: f32, y1: f32) -> Vec<[f32; 2]> { + fn rect(x0: f64, y0: f64, x1: f64, y1: f64) -> Vec<[f64; 2]> { vec![[x0, y0], [x1, y0], [x1, y1], [x0, y1]] } // Two nested rectangles, regardless of draw order, the resolution must be // deterministic and independent of which was drawn first. - fn nested(draw_order: bool) -> Vec> { + fn nested(draw_order: bool) -> Vec> { let big = rect(-10.0, -10.0, 10.0, 10.0); let small = rect(-5.0, -5.0, 5.0, 5.0); if draw_order { diff --git a/src/scene/boundary.rs b/src/scene/boundary.rs index 76e7d9b58..054ca4f38 100644 --- a/src/scene/boundary.rs +++ b/src/scene/boundary.rs @@ -3,7 +3,7 @@ use super::*; use std::collections::{HashMap, HashSet}; const NODE_EPS: f64 = 1.0e-6; -const PARAM_EPS: f64 = 1.0e-9; +const ORIENTATION_EPS: f64 = 1.0e-12; const AREA_EPS: f64 = 1.0e-10; #[derive(Clone, Copy, Debug)] @@ -43,6 +43,51 @@ impl Segment { } } +#[derive(Clone, Copy, Debug)] +struct SegmentAabb { + min_x: f64, + max_x: f64, + min_y: f64, + max_y: f64, +} + +impl SegmentAabb { + fn new(segment: Segment) -> Self { + Self { + min_x: segment.a.x.min(segment.b.x), + max_x: segment.a.x.max(segment.b.x), + min_y: segment.a.y.min(segment.b.y), + max_y: segment.a.y.max(segment.b.y), + } + } + + fn min(self, axis: SweepAxis) -> f64 { + match axis { + SweepAxis::X => self.min_x, + SweepAxis::Y => self.min_y, + } + } + + fn max(self, axis: SweepAxis) -> f64 { + match axis { + SweepAxis::X => self.max_x, + SweepAxis::Y => self.max_y, + } + } +} + +#[derive(Clone, Copy)] +enum SweepAxis { + X, + Y, +} + +enum SegmentIntersection { + None, + Point { a: f64, b: f64 }, + Overlap { a: [f64; 2], b: [f64; 2] }, +} + impl Scene { /// Build closed planar regions from the visible wire geometry. /// @@ -53,7 +98,7 @@ impl Scene { /// Curved entities participate through their already-tessellated WireModel /// geometry, so arcs, circles, ellipses and splines can take part in the /// boundary search without modifying the source entities. - pub fn hatch_boundary_outlines(&self) -> Vec> { + pub fn hatch_boundary_outlines(&self) -> Vec> { let wires = self.entity_wires(); let mut segments = Vec::::new(); @@ -87,7 +132,7 @@ impl Scene { build_planar_outlines(&segments) .into_iter() - .map(|ring| ring.into_iter().map(|p| [p.x as f32, p.y as f32]).collect()) + .map(|ring| ring.into_iter().map(|p| [p.x, p.y]).collect()) .collect() } } @@ -105,33 +150,57 @@ fn build_planar_outlines(segments: &[Segment]) -> Vec> { // Each original segment starts with its two endpoints as split positions. let mut cuts: Vec> = vec![vec![0.0, 1.0]; segments.len()]; - // Find all pairwise intersections. The AABB test avoids doing the more - // expensive intersection calculation for obviously unrelated segments. - for i in 0..segments.len() { - for j in (i + 1)..segments.len() { - if !segment_aabbs_overlap(segments[i], segments[j]) { + // Sweep the less-congested coordinate axis, then run the exact AABB and + // segment tests only for active candidates. Dense intersections remain + // output-sensitive, while spatially separated geometry avoids all-pairs + // work. + let aabbs: Vec = segments.iter().copied().map(SegmentAabb::new).collect(); + let axis = choose_sweep_axis(&aabbs); + let mut order: Vec = (0..segments.len()).collect(); + order.sort_by(|&a, &b| aabbs[a].min(axis).total_cmp(&aabbs[b].min(axis))); + let mut active = Vec::::new(); + + for i in order { + let current_min = aabbs[i].min(axis); + active.retain(|&j| aabbs[j].max(axis) + NODE_EPS >= current_min); + + for &j in &active { + if !segment_aabbs_overlap(aabbs[i], aabbs[j]) { continue; } - if let Some((ti, tj)) = segment_intersection(segments[i], segments[j]) { - cuts[i].push(ti.clamp(0.0, 1.0)); - cuts[j].push(tj.clamp(0.0, 1.0)); + match segment_intersection(segments[i], segments[j]) { + SegmentIntersection::None => {} + SegmentIntersection::Point { a, b } => { + cuts[i].push(a.clamp(0.0, 1.0)); + cuts[j].push(b.clamp(0.0, 1.0)); + } + SegmentIntersection::Overlap { a, b } => { + cuts[i].push(a[0].clamp(0.0, 1.0)); + cuts[i].push(a[1].clamp(0.0, 1.0)); + cuts[j].push(b[0].clamp(0.0, 1.0)); + cuts[j].push(b[1].clamp(0.0, 1.0)); + } } } + + active.push(i); } // Split every segment at all of its intersection parameters. let mut pieces = Vec::::new(); for (segment, params) in segments.iter().copied().zip(cuts.iter_mut()) { + let segment_len = segment.a.distance2(segment.b).sqrt(); + let param_merge_eps = (NODE_EPS / segment_len).min(1.0); params.sort_by(|a, b| a.total_cmp(b)); - params.dedup_by(|a, b| (*a - *b).abs() <= PARAM_EPS); + params.dedup_by(|a, b| (*a - *b).abs() <= param_merge_eps); for pair in params.windows(2) { let t0 = pair[0]; let t1 = pair[1]; - if t1 - t0 <= PARAM_EPS { + if (t1 - t0) * segment_len <= NODE_EPS { continue; } @@ -146,7 +215,7 @@ fn build_planar_outlines(segments: &[Segment]) -> Vec> { // Convert split segment endpoints to graph nodes. let mut nodes = Vec::::new(); - let mut node_map = HashMap::<(i64, i64), usize>::new(); + let mut node_map = HashMap::<(i64, i64), Vec>::new(); let mut edges = HashSet::<(usize, usize)>::new(); for piece in pieces { @@ -262,70 +331,142 @@ fn build_planar_outlines(segments: &[Segment]) -> Vec> { faces } -fn node_for_point(p: P2, nodes: &mut Vec, map: &mut HashMap<(i64, i64), usize>) -> usize { +fn node_for_point( + p: P2, + nodes: &mut Vec, + map: &mut HashMap<(i64, i64), Vec>, +) -> usize { let key = node_key(p); - if let Some(&index) = map.get(&key) { - return index; + for dx in -1_i64..=1 { + for dy in -1_i64..=1 { + let neighbor = (key.0.saturating_add(dx), key.1.saturating_add(dy)); + if let Some(indices) = map.get(&neighbor) { + if let Some(&index) = indices + .iter() + .find(|&&index| nodes[index].distance2(p) <= NODE_EPS * NODE_EPS) + { + return index; + } + } + } } let index = nodes.len(); nodes.push(p); - map.insert(key, index); + map.entry(key).or_default().push(index); index } fn node_key(p: P2) -> (i64, i64) { ( - (p.x / NODE_EPS).round() as i64, - (p.y / NODE_EPS).round() as i64, + (p.x / NODE_EPS).floor() as i64, + (p.y / NODE_EPS).floor() as i64, ) } -fn segment_aabbs_overlap(a: Segment, b: Segment) -> bool { - let a_min_x = a.a.x.min(a.b.x); - let a_max_x = a.a.x.max(a.b.x); - let a_min_y = a.a.y.min(a.b.y); - let a_max_y = a.a.y.max(a.b.y); +fn choose_sweep_axis(aabbs: &[SegmentAabb]) -> SweepAxis { + if normalized_interval_span(aabbs, SweepAxis::X) + <= normalized_interval_span(aabbs, SweepAxis::Y) + { + SweepAxis::X + } else { + SweepAxis::Y + } +} - let b_min_x = b.a.x.min(b.b.x); - let b_max_x = b.a.x.max(b.b.x); - let b_min_y = b.a.y.min(b.b.y); - let b_max_y = b.a.y.max(b.b.y); +fn normalized_interval_span(aabbs: &[SegmentAabb], axis: SweepAxis) -> f64 { + let min = aabbs + .iter() + .map(|aabb| aabb.min(axis)) + .fold(f64::INFINITY, f64::min); + let max = aabbs + .iter() + .map(|aabb| aabb.max(axis)) + .fold(f64::NEG_INFINITY, f64::max); + let extent = max - min; - a_max_x + NODE_EPS >= b_min_x - && b_max_x + NODE_EPS >= a_min_x - && a_max_y + NODE_EPS >= b_min_y - && b_max_y + NODE_EPS >= a_min_y + if !extent.is_finite() || extent <= NODE_EPS { + return f64::INFINITY; + } + + let total: f64 = aabbs + .iter() + .map(|aabb| aabb.max(axis) - aabb.min(axis) + NODE_EPS) + .sum(); + total / extent +} + +fn segment_aabbs_overlap(a: SegmentAabb, b: SegmentAabb) -> bool { + a.max_x + NODE_EPS >= b.min_x + && b.max_x + NODE_EPS >= a.min_x + && a.max_y + NODE_EPS >= b.min_y + && b.max_y + NODE_EPS >= a.min_y } /// Intersection parameters of two finite XY segments. -/// -/// Collinear overlaps intentionally return None. Shared endpoints and -/// T-junctions are handled by the normal non-parallel case and/or node merging. -fn segment_intersection(a: Segment, b: Segment) -> Option<(f64, f64)> { +fn segment_intersection(a: Segment, b: Segment) -> SegmentIntersection { let rx = a.b.x - a.a.x; let ry = a.b.y - a.a.y; let sx = b.b.x - b.a.x; let sy = b.b.y - b.a.y; - let cross = rx * sy - ry * sx; - - if cross.abs() <= 1.0e-12 { - return None; + let r2 = rx * rx + ry * ry; + let s2 = sx * sx + sy * sy; + if r2 <= f64::EPSILON || s2 <= f64::EPSILON { + return SegmentIntersection::None; } - + let r_len = r2.sqrt(); + let s_len = s2.sqrt(); + let cross = rx * sy - ry * sx; let qpx = b.a.x - a.a.x; let qpy = b.a.y - a.a.y; + if cross.abs() <= ORIENTATION_EPS * r_len * s_len { + if (qpx * ry - qpy * rx).abs() > NODE_EPS * r_len { + return SegmentIntersection::None; + } + + let b0_on_a = (qpx * rx + qpy * ry) / r2; + let b1_on_a = b0_on_a + (sx * rx + sy * ry) / r2; + let lo = b0_on_a.min(b1_on_a).max(0.0); + let hi = b0_on_a.max(b1_on_a).min(1.0); + let a_param_eps = NODE_EPS / r_len; + + if hi < lo - a_param_eps { + return SegmentIntersection::None; + } + + let lo = lo.clamp(0.0, 1.0); + let hi = hi.clamp(0.0, 1.0); + let lo_point = a.a.lerp(a.b, lo); + let hi_point = a.a.lerp(a.b, hi); + let b_lo = ((lo_point.x - b.a.x) * sx + (lo_point.y - b.a.y) * sy) / s2; + let b_hi = ((hi_point.x - b.a.x) * sx + (hi_point.y - b.a.y) * sy) / s2; + + if (hi - lo) * r_len <= NODE_EPS { + return SegmentIntersection::Point { + a: (lo + hi) * 0.5, + b: (b_lo + b_hi) * 0.5, + }; + } + + return SegmentIntersection::Overlap { + a: [lo, hi], + b: [b_lo, b_hi], + }; + } + let t = (qpx * sy - qpy * sx) / cross; let u = (qpx * ry - qpy * rx) / cross; + let t_eps = NODE_EPS / r_len; + let u_eps = NODE_EPS / s_len; - if t >= -PARAM_EPS && t <= 1.0 + PARAM_EPS && u >= -PARAM_EPS && u <= 1.0 + PARAM_EPS { - Some((t, u)) + if t >= -t_eps && t <= 1.0 + t_eps && u >= -u_eps && u <= 1.0 + u_eps { + SegmentIntersection::Point { a: t, b: u } } else { - None + SegmentIntersection::None } } @@ -334,12 +475,14 @@ fn signed_area(poly: &[P2]) -> f64 { return 0.0; } + let origin = poly[0]; let mut area = 0.0; - for i in 0..poly.len() { + for i in 1..poly.len() - 1 { let a = poly[i]; - let b = poly[(i + 1) % poly.len()]; - area += a.x * b.y - b.x * a.y; + let b = poly[i + 1]; + area += (a.x - origin.x) * (b.y - origin.y) + - (b.x - origin.x) * (a.y - origin.y); } area * 0.5 @@ -395,4 +538,87 @@ mod tests { assert!(faces.is_empty()); } + + #[test] + fn collinear_overlap_closes_rectangle() { + let segments = vec![ + s(0.0, 0.0, 7.0, 0.0), + s(3.0, 0.0, 10.0, 0.0), + s(10.0, 0.0, 10.0, 10.0), + s(10.0, 10.0, 0.0, 10.0), + s(0.0, 10.0, 0.0, 0.0), + ]; + + let faces = build_planar_outlines(&segments); + + assert_eq!(faces.len(), 1); + assert!((signed_area(&faces[0]) - 100.0).abs() < 1.0e-6); + } + + #[test] + fn node_merge_crosses_bucket_boundary() { + let mut nodes = Vec::new(); + let mut map = HashMap::new(); + let a = node_for_point(P2::new(NODE_EPS * 0.99, 0.0), &mut nodes, &mut map); + let b = node_for_point(P2::new(NODE_EPS * 1.01, 0.0), &mut nodes, &mut map); + + assert_eq!(a, b); + } + + #[test] + fn node_merge_rejects_distant_diagonal_points() { + let mut nodes = Vec::new(); + let mut map = HashMap::new(); + let a = node_for_point( + P2::new(NODE_EPS * 0.51, NODE_EPS * 0.51), + &mut nodes, + &mut map, + ); + let b = node_for_point( + P2::new(NODE_EPS * 1.49, NODE_EPS * 1.49), + &mut nodes, + &mut map, + ); + + assert_ne!(a, b); + } + + #[test] + fn signed_area_is_stable_at_large_coordinates() { + let base = 1.0e12; + let poly = vec![ + P2::new(base, base), + P2::new(base + 3.0, base), + P2::new(base + 3.0, base + 4.0), + P2::new(base, base + 4.0), + ]; + + assert!((signed_area(&poly) - 12.0).abs() < 1.0e-9); + } + + #[test] + fn long_segments_keep_world_scale_cuts() { + let segments = vec![ + s(0.0, 0.0, 1.0e12, 0.0), + s(0.0, 1.0, 1.0e12, 1.0), + s(1.0, -1.0, 1.0, 2.0), + s(2.0, -1.0, 2.0, 2.0), + ]; + + let faces = build_planar_outlines(&segments); + + assert_eq!(faces.len(), 1); + assert!((signed_area(&faces[0]) - 1.0).abs() < 1.0e-9); + } + + #[test] + fn long_near_collinear_segments_overlap_within_node_tolerance() { + let a = s(0.0, 0.0, 1.0e9, 1.0e-3); + let b = s(5.0e8, 5.0e-4 + 5.0e-7, 1.5e9, 1.5e-3 + 5.0e-7); + + assert!(matches!( + segment_intersection(a, b), + SegmentIntersection::Overlap { .. } + )); + } } diff --git a/src/scene/mod.rs b/src/scene/mod.rs index 90403b379..52acf6d89 100644 --- a/src/scene/mod.rs +++ b/src/scene/mod.rs @@ -4590,11 +4590,27 @@ impl Scene { } /// Collect closed polygon outlines (world XY) from the current layout. - pub fn closed_outlines(&self) -> Vec> { + pub fn closed_outlines(&self) -> Vec> { self.entity_wires() .iter() .filter_map(|wire| { - let pts = &wire.points; + let pts: Vec<[f64; 2]> = wire + .points + .iter() + .copied() + .enumerate() + .filter_map(|(index, high)| { + if !high[0].is_finite() || !high[1].is_finite() { + return None; + } + let low = wire.points_low.get(index).copied().unwrap_or([0.0; 3]); + let point = [ + high[0] as f64 + low[0] as f64, + high[1] as f64 + low[1] as f64, + ]; + point.iter().all(|value| value.is_finite()).then_some(point) + }) + .collect(); if pts.len() < 4 { return None; } @@ -4612,12 +4628,8 @@ impl Scene { // previous one, so consumers (point-in-polygon, the hatch / // boundary commands) see one vertex per corner — not the doubled // ring that otherwise shows two grips at every corner. - let mut ring: Vec<[f32; 2]> = Vec::with_capacity(pts.len()); - for p in pts { - if !p[0].is_finite() || !p[1].is_finite() { - continue; - } - let q = [p[0], p[1]]; + let mut ring: Vec<[f64; 2]> = Vec::with_capacity(pts.len()); + for q in pts { if let Some(&last) = ring.last() { if (last[0] - q[0]).abs() < 1e-4 && (last[1] - q[1]).abs() < 1e-4 { continue;