diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index 1b050960..d6668d5c 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -31,3 +31,9 @@ Generated-in-code assets are original and have no file to list: | crates/mcp/src/server.rs, crates/mcp/src/backend.rs (JSON-RPC framing, MCP lifecycle, control-port bridge) | CADCraft `crates/mcp/src/{server,backend}.rs` | MIT OR Apache-2.0 | | crates/mcp/src/tools.rs (tool schema helpers and argument checking) | CADCraft `crates/mcp/src/tools.rs`, GridCraft `crates/mcp/src/tools.rs` | MIT OR Apache-2.0 | | LICENSE-MIT, LICENSE-APACHE (licence text) | CADCraft | — | + +## Vendored third-party code + +| Directory | Source | Licence | Changes | +|---|---|---|---| +| vendor/truck-shapeops | truck-shapeops 0.4.0, https://github.com/ricosjp/truck (RICOS Co. Ltd.) | Apache-2.0 | robust meshing and nearest-parameter fallback in booleans; trace macro (vendor/README.md) | diff --git a/Cargo.lock b/Cargo.lock index 24f3617e..cf0f8424 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3502,8 +3502,6 @@ dependencies = [ [[package]] name = "truck-shapeops" version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "fd6e83ea4d36d00229d38ab985a991e32f58c3ed2971083011ac66fb57e26d5c" dependencies = [ "derive_more", "itertools 0.13.0", diff --git a/Cargo.toml b/Cargo.toml index 339c4fa0..326b2da3 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,8 @@ [workspace] resolver = "3" members = ["crates/*", "apps/*", "xtask"] +# Third-party code kept with local fixes (see vendor/README.md); not linted as ours. +exclude = ["vendor"] [workspace.package] version = "0.1.0" @@ -68,3 +70,7 @@ opt-level = 3 lto = "thin" codegen-units = 16 debug = 0 + +# truck-shapeops with SolveCraft fixes (Apache-2.0, see vendor/truck-shapeops and ATTRIBUTION.md). +[patch.crates-io] +truck-shapeops = { path = "vendor/truck-shapeops" } diff --git a/NOTICE b/NOTICE index 27f1cf20..36bf1a5c 100644 --- a/NOTICE +++ b/NOTICE @@ -8,7 +8,9 @@ Third-party material Software dependencies are listed in Cargo.lock. The geometry kernel builds on the truck crates (truck-modeling, truck-shapeops, truck-meshalgo, truck-stepio, truck-topology and their dependencies, https://github.com/ricosjp/truck), licensed under the Apache License 2.0 and used -unmodified as libraries. +as libraries. truck-shapeops 0.4.0 is vendored in vendor/truck-shapeops with SolveCraft +modifications (listed in vendor/README.md and marked in the source); its licence is in +vendor/truck-shapeops/LICENSE. Copyright (c) RICOS Co. Ltd. Some code is adapted from sibling projects in the storytold family (CADCraft, GridCraft), which are licensed MIT OR Apache-2.0; ATTRIBUTION.md lists the files. diff --git a/ROADMAP.md b/ROADMAP.md index 26879bf1..72032b78 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -9,7 +9,7 @@ features → parametric timeline) as an open, pure-Rust, agent-drivable applicat |---|---| | Where are we? | **M0 done; M1/M2 in progress.** Sketch with solver, extrude (taper, two-sided, through all), revolve, fillet/chamfer, holes, patterns, mirror, shell, draft, loft, sweep, split, construction planes, parameters that rebuild the timeline, STEP/STL export, desktop and browser (wasm) app with Fusion-style selection (hover/selected highlights, origin planes, click-to-pick sketch planes, multi-select, window/crossing box selection, animated view cube), CLI, MCP server, oracle harness. | | Command parity (SOLID + SKETCH toolbar) | **57 / 291 (20%)** — [docs/parity.md](docs/parity.md) | -| Fusion oracle | **23 / 29 parts match** Fusion's measurements — [docs/oracle.md](docs/oracle.md) | +| Fusion oracle | **25 / 29 parts match** Fusion's measurements — [docs/oracle.md](docs/oracle.md) | | Tests | 85 (solver, profiles, kernel booleans/blends/measures, expressions, timeline, file formats, camera, engine end-to-end, hostile-input fuzz over every command) | | Gates | `cargo xtask ci`: fmt, clippy -D warnings, tests, asset attribution, layering, wasm32 build — green | | Weighted parity estimate | **≈ 5%** of Fusion's Design workspace by importance (sketch + basic solids are the core, but surfaces, assemblies, sheet metal, CAM, drawings are untouched) | @@ -37,8 +37,9 @@ features → parametric timeline) as an open, pure-Rust, agent-drivable applicat - Fillets and chamfers: straight convex edges between planar faces with perpendicular planar end faces; two blended edges may not share a corner yet. -- Booleans of bodies with coincident faces fail unless they come from extrudes (where the tool is - extended automatically); fully internal voids are not supported. +- Booleans of curved bodies with coincident faces can fail (planar ones have an exact fallback); + fully internal voids are not supported. Spheres are built from six pole-free patches so they + combine reliably. - Mass properties come from fine tessellation (curved faces within ~1e-4 relative). - STEP import reads solids (B-rep with analytic and B-spline geometry), units, product names, colours and assemblies; offset surfaces and pcurve-only edges are not read yet, and assembly diff --git a/crates/kernel/src/body.rs b/crates/kernel/src/body.rs index 9cd6e5ea..6f9320f9 100644 --- a/crates/kernel/src/body.rs +++ b/crates/kernel/src/body.rs @@ -37,17 +37,16 @@ pub struct FaceInfo { type MeshedShell = truck_topology::Shell, Option>; /// Triangulate each shell (not the solid: truck's solid constructor panics on the open or -/// non-manifold shells an import can produce). Faces whose boundary is not exactly on the -/// surface (imported geometry within the file's tolerance) fail the exact parameter search; -/// such shells are meshed again projecting onto the surface. +/// non-manifold shells an import can produce). Boundaries that lie on their surface only within +/// tolerance (imports, boolean results) are projected onto it. fn mesh_shells(solid: &Solid, tol: f64) -> Vec { solid .boundaries() .iter() - .map(|sh| { - let m = sh.triangulation(tol); - if m.face_iter().any(|f| f.surface().is_none()) { sh.robust_triangulation(tol) } else { m } - }) + // The robust search tries the exact one first, so exact geometry costs nothing extra; a + // face can also come back as an empty mesh rather than none (boolean results with + // approximate intersection curves), so always use it. + .map(|sh| sh.robust_triangulation(tol)) .collect() } diff --git a/crates/kernel/src/build.rs b/crates/kernel/src/build.rs index 224956fe..bfa4018d 100644 --- a/crates/kernel/src/build.rs +++ b/crates/kernel/src/build.rs @@ -140,6 +140,14 @@ pub fn revolve(plane: &Plane, regions: &[Region2], axis_origin: Vec2, axis_dir: .collect(), }; let r = &Region2 { outer: snap_loop(&r.outer), holes: r.holes.iter().map(snap_loop).collect() }; + // A half disc on the axis revolved fully is a sphere: build it without poles. + if angle.abs() >= std::f64::consts::TAU - 1e-9 + && r.holes.is_empty() + && let Some((c, rad)) = half_disc(&r.outer, axis_origin, d) + { + out.push(sphere_patches(plane.to_world(c), rad)?); + continue; + } // A full revolution of a profile with one edge on the axis: sweep the rest of the loop // as a cone-like shell (no zero-area faces, which booleans can't handle). if angle.abs() >= std::f64::consts::TAU - 1e-9 @@ -191,19 +199,7 @@ pub fn cylinder(base: Vec3, axis: Vec3, radius: f64, height: f64) -> Result Result { - let radius = finite(radius, "radius")?; - if radius <= 1e-6 { - return Err(KernelError::Invalid("radius must be positive".into())); - } - let plane = Plane { origin: center, ..Plane::XZ }; - let half = Loop2 { - segs: vec![ - Seg2::Arc { center: Vec2::ZERO, radius, start: -std::f64::consts::FRAC_PI_2, sweep: std::f64::consts::PI }, - Seg2::Line { a: Vec2::new(0.0, radius), b: Vec2::new(0.0, -radius) }, - ], - }; - let mut v = revolve(&plane, &[Region2 { outer: half, holes: vec![] }], Vec2::ZERO, Vec2::Y, std::f64::consts::TAU)?; - v.pop().ok_or_else(|| KernelError::Failed("sphere".into())) + sphere_patches(center, radius) } /// Torus about the Z axis through `center`. @@ -313,6 +309,23 @@ pub fn extrude_tapered(plane: &Plane, region: &Region2, length: f64, dir_sign: f } /// Revolve 360° a loop that has exactly one line segment on the axis, without degenerate faces. +/// A loop made of a half circle and its diameter on the axis: (centre, radius). +fn half_disc(lp: &Loop2, axis_origin: Vec2, d: Vec2) -> Option<(Vec2, f64)> { + let on_axis = |p: Vec2| d.cross(p - axis_origin).abs() < 1e-7 * (1.0 + p.len()); + let mut arcs = Vec::new(); + for s in &lp.segs { + match *s { + Seg2::Line { a, b } if on_axis(a) && on_axis(b) => {} + Seg2::Arc { center, radius, sweep, .. } if on_axis(center) => arcs.push((center, radius, sweep)), + _ => return None, + } + } + let (c, r, _) = *arcs.first()?; + let total: f64 = arcs.iter().map(|a| a.2.abs()).sum(); + let same = arcs.iter().all(|(c2, r2, _)| c2.dist(c) < 1e-7 * (1.0 + r) && (r2 - r).abs() < 1e-7 * (1.0 + r)); + (same && (total - std::f64::consts::PI).abs() < 1e-6).then_some((c, r)) +} + fn revolve_touching_axis(plane: &Plane, lp: &Loop2, axis_origin: Vec2, d: Vec2) -> Result> { let on_axis = |p: Vec2| d.cross(p - axis_origin).abs() < 1e-9; let n = lp.segs.len(); @@ -654,3 +667,77 @@ fn uniform_wire(plane: &Plane, lp: &Loop2) -> Result { } Ok(edges.into()) } + +/// A sphere as six patches (the cube's faces projected onto it) on revolved surfaces whose poles +/// and seams lie outside each patch. No face has a singular point, which the boolean needs; the +/// patches are one analytic sphere, so measured topology counts it as one face. +pub fn sphere_patches(center: Vec3, radius: f64) -> Result { + sphere_patches_rotated(center, radius, [Vec3::X, Vec3::Y, Vec3::Z]) +} + +/// [`sphere_patches`] with the patch layout turned to the frame `rot` (orthonormal columns). +pub fn sphere_patches_rotated(center: Vec3, radius: f64, rot: [Vec3; 3]) -> Result { + let radius = finite(radius, "radius")?; + let turn = |v: Vec3| rot[0] * v.x + rot[1] * v.y + rot[2] * v.z; + if radius <= 1e-6 || !center.is_finite() { + return Err(KernelError::Invalid("radius must be positive".into())); + } + let k = radius / 3f64.sqrt(); + let corner = |sx: f64, sy: f64, sz: f64| center + turn(Vec3::new(sx * k, sy * k, sz * k)); + // Cube faces: outward direction, rotation axis for the patch's surface (perpendicular to it), + // and the four corner sign triples counter-clockwise seen from outside. + let faces: [(Vec3, Vec3, [[f64; 3]; 4]); 6] = [ + (Vec3::X, Vec3::Z, [[1.0, -1.0, -1.0], [1.0, 1.0, -1.0], [1.0, 1.0, 1.0], [1.0, -1.0, 1.0]]), + (-Vec3::X, Vec3::Z, [[-1.0, 1.0, -1.0], [-1.0, -1.0, -1.0], [-1.0, -1.0, 1.0], [-1.0, 1.0, 1.0]]), + (Vec3::Y, Vec3::Z, [[1.0, 1.0, -1.0], [-1.0, 1.0, -1.0], [-1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]), + (-Vec3::Y, Vec3::Z, [[-1.0, -1.0, -1.0], [1.0, -1.0, -1.0], [1.0, -1.0, 1.0], [-1.0, -1.0, 1.0]]), + (Vec3::Z, Vec3::X, [[-1.0, -1.0, 1.0], [1.0, -1.0, 1.0], [1.0, 1.0, 1.0], [-1.0, 1.0, 1.0]]), + (-Vec3::Z, Vec3::X, [[-1.0, 1.0, -1.0], [1.0, 1.0, -1.0], [1.0, -1.0, -1.0], [-1.0, -1.0, -1.0]]), + ]; + guard("sphere", || { + let key = |s: [f64; 3]| ((s[0] > 0.0) as usize) | (((s[1] > 0.0) as usize) << 1) | (((s[2] > 0.0) as usize) << 2); + let mut verts: std::collections::HashMap = std::collections::HashMap::new(); + let mut vert = |s: [f64; 3]| verts.entry(key(s)).or_insert_with(|| builder::vertex(p3(corner(s[0], s[1], s[2])))).clone(); + let mut edges: std::collections::HashMap<(usize, usize), mt::Edge> = std::collections::HashMap::new(); + let mut out: Vec = Vec::new(); + for (dir, axis, cs) in faces { + let (dir, axis) = (turn(dir), turn(axis)); + let mut wire: Vec = Vec::new(); + for i in 0..4 { + let (a, b) = (cs[i], cs[(i + 1) % 4]); + let (ka, kb) = (key(a), key(b)); + let e = match edges.get(&(ka.min(kb), ka.max(kb))) { + Some(e) => e.clone(), + None => { + let (lo, hi) = if ka < kb { (a, b) } else { (b, a) }; + let (va, vb) = (vert(lo), vert(hi)); + let mid = (corner(lo[0], lo[1], lo[2]) + corner(hi[0], hi[1], hi[2])) * 0.5 - center; + let transit = center + mid.normalized().unwrap_or(dir) * radius; + let e = builder::circle_arc(&va, &vb, p3(transit)); + edges.insert((ka.min(kb), ka.max(kb)), e.clone()); + e + } + }; + wire.push(if ka < kb { e } else { e.inverse() }); + } + // The meridian runs pole to pole through the side away from the patch. + let (np, sp, back) = (center + axis * radius, center - axis * radius, center - dir * radius); + let meridian = builder::circle_arc(&builder::vertex(p3(np)), &builder::vertex(p3(sp)), p3(back)); + let mut surface = + mt::Surface::RevolutedCurve(mt::Processor::new(mt::RevolutedCurve::by_revolution(meridian.oriented_curve(), p3(center), v3(axis)))); + // Outward normal at the patch centre. + let probe = p3(center + dir * radius); + use mt::{ParametricSurface3D, SearchParameter}; + if let Some((u, v)) = surface.search_parameter(probe, mt::SPHint2D::None, 100) { + let n = surface.normal(u, v); + if mt::InnerSpace::dot(n, v3(dir)) < 0.0 { + mt::Invertible::invert(&mut surface); + } + } + out.push(mt::Face::new(vec![wire.into()], surface)); + } + let shell: mt::Shell = out.into(); + let solid = Solid::try_new(vec![shell]).map_err(|e| KernelError::Failed(format!("sphere: {e}")))?; + Body::new(solid) + }) +} diff --git a/crates/kernel/src/tests.rs b/crates/kernel/src/tests.rs index 9a09a279..dfbd890d 100644 --- a/crates/kernel/src/tests.rs +++ b/crates/kernel/src/tests.rs @@ -246,6 +246,20 @@ fn planar_booleans_with_coincident_faces() { assert!(rel(measure(&i).unwrap().volume, 250.0) < 1e-9); } +#[test] +fn sphere_booleans() { + let s = sphere(Vec3::ZERO, 15.0).unwrap(); + let m = measure(&s).unwrap(); + assert!(rel(m.volume, 4.0 / 3.0 * PI * 3375.0) < 1e-3, "{}", m.volume); + assert_eq!((m.merged.faces, m.merged.edges), (1, 0)); + // A ball cut out of the top face of a slab (oracle 14). + let b = box_solid(Vec3::new(-25.0, -25.0, -20.0), Vec3::new(25.0, 25.0, 0.0)).unwrap(); + let c = boolean(&b, &s, BoolOp::Cut).unwrap().unwrap(); + let m = measure(&c).unwrap(); + assert!(rel(m.volume, 50.0 * 50.0 * 20.0 - 2.0 / 3.0 * PI * 3375.0) < 1e-3, "{}", m.volume); + assert_eq!((m.merged.faces, m.merged.edges, m.merged.vertices), (7, 13, 9)); +} + #[test] #[ignore] fn debug_convergence() { @@ -456,3 +470,132 @@ fn debug_cbore_seq() { println!(" swapped: {:?}", boolean(&c2, &t1, BoolOp::Cut).map(|o| o.map(|x| measure(&x).unwrap().volume))); } } + +#[test] +#[ignore] +fn debug_truck_sphere() { + use truck_modeling::{self as mt, builder}; + let mk = |c: Vec3, r: f64, axis: Vec3| -> Body { + let ax = axis.normalized().unwrap(); + let side = ax.any_perp(); + let v0 = builder::vertex(crate::body::p3(c + ax * r)); + let arc = builder::rsweep(&v0, crate::body::p3(c), crate::body::v3(side), mt::Rad(PI)); + let shell = builder::rsweep(&arc, crate::body::p3(c), crate::body::v3(ax), mt::Rad(2.0 * PI)); + Body::new(mt::Solid::new(vec![shell])).unwrap() + }; + let b = box_solid(Vec3::new(-25.0, -25.0, -20.0), Vec3::new(25.0, 25.0, 0.0)).unwrap(); + for (name, axis) in [("z", Vec3::Z), ("x", Vec3::X), ("diag", Vec3::new(0.3, 0.5, 0.8))] { + let sp = mk(Vec3::ZERO, 15.0, axis); + println!("{name}: faces {} vol {}", sp.face_count(), measure(&sp).unwrap().volume); + for op in [BoolOp::Cut, BoolOp::Intersect, BoolOp::Union] { + let r = guard("t", || { + let (sa, mut sb) = (b.deep_copy(), sp.deep_copy()); + let res = match op { + BoolOp::Cut => { + sb.not(); + truck_shapeops::and(&sa, &sb, 0.01) + } + BoolOp::Intersect => truck_shapeops::and(&sa, &sb, 0.01), + BoolOp::Union => truck_shapeops::or(&sa, &sb, 0.01), + }; + Ok(res.map(|s| s.face_iter().count())) + }); + println!(" {op:?}: {r:?}"); + } + } +} + +#[test] +#[ignore] +fn debug_sphere_patches() { + let sp = crate::build::sphere_patches(Vec3::ZERO, 15.0).unwrap(); + let m = measure(&sp).unwrap(); + println!("vol {} want {} area {} merged {:?}", m.volume, 4.0 / 3.0 * PI * 3375.0, m.area, m.merged); + let b = box_solid(Vec3::new(-25.0, -25.0, -20.0), Vec3::new(25.0, 25.0, 0.0)).unwrap(); + for op in [BoolOp::Cut, BoolOp::Intersect, BoolOp::Union] { + let r = boolean(&b, &sp, op).map(|o| { + o.map(|x| { + let m = measure(&x).unwrap(); + (m.volume, m.merged) + }) + }); + println!("{op:?}: {r:?}"); + } + let c = box_solid(Vec3::new(-20.0, -20.0, -20.0), Vec3::new(20.0, 20.0, 20.0)).unwrap(); + let a = 0.3f64; + let (ca, sa) = (a.cos(), a.sin()); + let b2 = 0.5f64; + let (cb, sb) = (b2.cos(), b2.sin()); + // Rz(a) * Rx(b) + let rot = [Vec3::new(ca, sa, 0.0), Vec3::new(-sa * cb, ca * cb, sb), Vec3::new(sa * sb, -ca * sb, cb)]; + for (name, s2) in [ + ("plain", crate::build::sphere_patches(Vec3::ZERO, 25.0).unwrap()), + ("turned", crate::build::sphere_patches_rotated(Vec3::ZERO, 25.0, rot).unwrap()), + ] { + println!("{name} vol {}", measure(&s2).unwrap().volume); + println!( + "{name} cube and sphere: {:?}", + boolean(&c, &s2, BoolOp::Intersect).map(|o| o.map(|x| { + let m = measure(&x).unwrap(); + (m.volume, m.merged) + })) + ); + } +} + +#[test] +#[ignore] +fn debug_sphere_cases() { + let s = crate::build::sphere_patches(Vec3::ZERO, 25.0).unwrap(); + let cases = [ + ("one plane", Vec3::new(-100.0, -100.0, -100.0), Vec3::new(100.0, 100.0, 20.0)), + ("two planes", Vec3::new(-100.0, -100.0, -20.0), Vec3::new(100.0, 100.0, 20.0)), + ("three", Vec3::new(-100.0, -100.0, -20.0), Vec3::new(20.0, 100.0, 20.0)), + ("four", Vec3::new(-100.0, -20.0, -20.0), Vec3::new(20.0, 20.0, 20.0)), + ("five", Vec3::new(-20.0, -20.0, -20.0), Vec3::new(20.0, 20.0, 100.0)), + ("cube", Vec3::new(-20.0, -20.0, -20.0), Vec3::new(20.0, 20.0, 20.0)), + ("cube big sphere-ish", Vec3::new(-24.0, -24.0, -24.0), Vec3::new(24.0, 24.0, 24.0)), + ("cube off", Vec3::new(-20.3, -19.7, -20.1), Vec3::new(19.9, 20.2, 19.8)), + ]; + for (name, lo, hi) in cases { + let b = box_solid(lo, hi).unwrap(); + for op in [BoolOp::Intersect, BoolOp::Cut] { + let r = boolean(&b, &s, op).map(|o| o.map(|x| measure(&x).unwrap().volume)); + println!("{name} {op:?}: {:?}", r.map_err(|e| e.to_string().chars().take(60).collect::())); + } + } +} + +#[test] +#[ignore] +fn debug_sphere_rotations() { + let b = box_solid(Vec3::new(-20.0, -20.0, -20.0), Vec3::new(20.0, 20.0, 20.0)).unwrap(); + for k in 0..8 { + let (a, c) = (0.37 * k as f64 + 0.1, 0.61 * k as f64 + 0.2); + let (ca, sa, cc, sc) = (a.cos(), a.sin(), c.cos(), c.sin()); + let rot = [Vec3::new(ca, sa, 0.0), Vec3::new(-sa * cc, ca * cc, sc), Vec3::new(sa * sc, -ca * sc, cc)]; + let s = crate::build::sphere_patches_rotated(Vec3::ZERO, 25.0, rot).unwrap(); + let r = guard("t", || { + let res = truck_shapeops::and(&b.deep_copy(), &s.deep_copy(), 0.01); + Ok(res.map(|x| x.face_iter().count())) + }); + println!("rot {k}: {r:?}"); + } +} + +#[test] +#[ignore] +fn debug_cube_sphere_measure() { + let b = box_solid(Vec3::new(-20.0, -20.0, -20.0), Vec3::new(20.0, 20.0, 20.0)).unwrap(); + let s = crate::build::sphere_patches(Vec3::ZERO, 25.0).unwrap(); + let r = guard("t", || Ok(truck_shapeops::and(&b.deep_copy(), &s.deep_copy(), 0.05))).unwrap().unwrap(); + let body = Body::new(r).unwrap(); + for tol in [0.5, 0.1, 0.02] { + let m = body.tessellate(tol).unwrap(); + let mm = m.measure(); + println!("tol {tol}: vol {} area {} tris {}", mm.volume, mm.area, m.triangles.len()); + } + for f in body.faces(0.05).unwrap() { + println!("face {} area {:.2} c {:?} planar {:?}", f.index, f.area, f.centroid, f.plane_normal.is_some()); + } +} diff --git a/crates/kernel/src/topo.rs b/crates/kernel/src/topo.rs index f98389e9..6834bc0e 100644 --- a/crates/kernel/src/topo.rs +++ b/crates/kernel/src/topo.rs @@ -135,23 +135,10 @@ fn classify(pts: &[Vec3], nrm: &[Vec3], tol: f64, id: usize) -> Surf { } } } - // Sphere: points equidistant from the point where normals meet. - let c_est = { - let mut acc = Vec3::ZERO; - let mut k = 0.0; - for (p, n) in pts.iter().zip(nrm) { - // Assume radius from the first pair. - let r0 = (p0 - *p).len() / ((n0 - *n).len().max(1e-12)); - if r0.is_finite() && r0 > 0.0 && (n0 - *n).len() > 0.1 { - acc += *p - *n * r0; - k += 1.0; - } - } - (k > 0.0).then(|| acc / k) - }; - if let Some(c) = c_est { - let r = pts.iter().map(|p| p.dist(c)).sum::() / pts.len().max(1) as f64; - if pts.iter().all(|p| (p.dist(c) - r).abs() < tol.max(r * 2e-3)) { + // Sphere: algebraic least-squares fit |p|² = 2c·p + k, then every point at the same distance. + if let Some((c, r)) = fit_sphere(pts) { + let spread = pts.iter().map(|p| p.dist(p0)).fold(0.0, f64::max); + if spread > tol && pts.iter().all(|p| (p.dist(c) - r).abs() < tol.max(r * 2e-3)) && r < spread * 1e4 { return Surf::Sphere { c, r }; } } @@ -278,6 +265,47 @@ pub fn seam_flags(b: &Body, mesh: &Mesh) -> Vec { .collect() } +/// Least-squares sphere through points: centre and radius. +fn fit_sphere(pts: &[Vec3]) -> Option<(Vec3, f64)> { + if pts.len() < 4 { + return None; + } + // Centre the data for conditioning. + let m = pts.iter().fold(Vec3::ZERO, |a, p| a + *p) / pts.len() as f64; + let mut a = [[0.0f64; 5]; 4]; + for p in pts { + let q = *p - m; + let row = [2.0 * q.x, 2.0 * q.y, 2.0 * q.z, 1.0]; + let rhs = q.dot(q); + for i in 0..4 { + for j in 0..4 { + a[i][j] += row[i] * row[j]; + } + a[i][4] += row[i] * rhs; + } + } + for col in 0..4 { + let piv = (col..4).max_by(|x, y| a[*x][col].abs().total_cmp(&a[*y][col].abs()))?; + if a[piv][col].abs() < 1e-18 { + return None; + } + a.swap(col, piv); + for r in 0..4 { + if r != col { + let f = a[r][col] / a[col][col]; + let prow = a[col]; + for (x, v) in a[r].iter_mut().zip(prow) { + *x -= f * v; + } + } + } + } + let c = Vec3::new(a[0][4] / a[0][0], a[1][4] / a[1][1], a[2][4] / a[2][2]); + let k = a[3][4] / a[3][3]; + let r2 = k + c.dot(c); + (r2 > 0.0 && r2.is_finite()).then(|| (c + m, r2.sqrt())) +} + fn classify_faces(mesh: &Mesh, nf: usize, tol: f64) -> Vec { let mut pts: Vec> = vec![Vec::new(); nf]; let mut nrm: Vec> = vec![Vec::new(); nf]; diff --git a/vendor/README.md b/vendor/README.md new file mode 100644 index 00000000..ff095c57 --- /dev/null +++ b/vendor/README.md @@ -0,0 +1,15 @@ +# Vendored third-party crates + +Copies of third-party crates that SolveCraft builds with local fixes, wired in with +`[patch.crates-io]` in the workspace `Cargo.toml`. They keep their own licence; changed lines +are marked `SolveCraft:` in the source. + +## truck-shapeops 0.4.0 (Apache-2.0, © RICOS Co. Ltd., https://github.com/ricosjp/truck) + +Changes: +- `transversal/integrate/mod.rs`: mesh both shells with `robust_triangulation` (faces made by + earlier booleans have boundaries that lie on their surface only within tolerance). +- `transversal/divide_face/mod.rs`: when projecting a loop onto a face's surface, fall back to + the nearest parameter when the exact search fails (intersection curves are approximations). +- `lib.rs`: `tr!` reports where an operation gives up when `SHAPEOPS_TRACE` is set; compiler + warnings no longer fail the build. diff --git a/vendor/truck-shapeops/.cargo-ok b/vendor/truck-shapeops/.cargo-ok new file mode 100644 index 00000000..5f8b7958 --- /dev/null +++ b/vendor/truck-shapeops/.cargo-ok @@ -0,0 +1 @@ +{"v":1} \ No newline at end of file diff --git a/vendor/truck-shapeops/Cargo.toml b/vendor/truck-shapeops/Cargo.toml new file mode 100644 index 00000000..8a20ea1f --- /dev/null +++ b/vendor/truck-shapeops/Cargo.toml @@ -0,0 +1,78 @@ +# THIS FILE IS AUTOMATICALLY GENERATED BY CARGO +# +# When uploading crates to the registry Cargo will automatically +# "normalize" Cargo.toml files for maximal compatibility +# with all versions of Cargo and also rewrite `path` dependencies +# to registry (e.g., crates.io) dependencies. +# +# If you are reading this file be aware that the original Cargo.toml +# will likely look very different (and much more reasonable). +# See Cargo.toml.orig for the original contents. + +[package] +edition = "2021" +name = "truck-shapeops" +version = "0.4.0" +build = false +autobins = false +autoexamples = false +autotests = false +autobenches = false +description = "Provides boolean operations to Solid" +homepage = "https://github.com/ricosjp/truck" +readme = "README.md" +license = "Apache-2.0" +repository = "https://github.com/ricosjp/truck" + +[lib] +name = "truck_shapeops" +path = "src/lib.rs" + +[[example]] +name = "punched-cube-shapeops" +path = "examples/punched-cube-shapeops.rs" + +[dependencies.derive_more] +version = "0.99.18" + +[dependencies.itertools] +version = "0.13.0" + +[dependencies.rustc-hash] +version = "2.0.0" + +[dependencies.truck-base] +version = "0.5.0" + +[dependencies.truck-geometry] +version = "0.5.0" + +[dependencies.truck-geotrait] +version = "0.4.0" + +[dependencies.truck-meshalgo] +version = "0.4.0" + +[dependencies.truck-stepio] +version = "0.3.0" +optional = true + +[dependencies.truck-topology] +version = "0.6.0" + +[dev-dependencies.rand] +version = "0.8.5" + +[dev-dependencies.serde_json] +version = "1.0.128" + +[dev-dependencies.truck-modeling] +version = "0.6.0" + +[features] +default = ["step-test"] +step-test = ["truck-stepio"] + +[target.'cfg(target_arch = "wasm32")'.dev-dependencies.getrandom] +version = "0.2.15" +features = ["js"] diff --git a/vendor/truck-shapeops/LICENSE b/vendor/truck-shapeops/LICENSE new file mode 100644 index 00000000..f433b1a5 --- /dev/null +++ b/vendor/truck-shapeops/LICENSE @@ -0,0 +1,177 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS diff --git a/vendor/truck-shapeops/README.md b/vendor/truck-shapeops/README.md new file mode 100644 index 00000000..5937526a --- /dev/null +++ b/vendor/truck-shapeops/README.md @@ -0,0 +1,11 @@ +# truck-shapeops + +[![Crates.io](https://img.shields.io/crates/v/truck-shapeops.svg)](https://crates.io/crates/truck-shapeops) [![Docs.rs](https://docs.rs/truck-shapeops/badge.svg)](https://docs.rs/truck-shapeops) + +Crate for operationg shapes. Provides boolean operations to Solid, and shape healing for importing shapes from other CAD systems. + +## Sample Codes + +### punched-cube-shapeops + +A punched cube generated by shape operator `AND`. diff --git a/vendor/truck-shapeops/examples/punched-cube-shapeops.rs b/vendor/truck-shapeops/examples/punched-cube-shapeops.rs new file mode 100644 index 00000000..4e5fc759 --- /dev/null +++ b/vendor/truck-shapeops/examples/punched-cube-shapeops.rs @@ -0,0 +1,30 @@ +//! A punched cube generated by shape operator `AND`. + +use truck_modeling::*; + +fn main() { + let v = builder::vertex(Point3::origin()); + let e = builder::tsweep(&v, Vector3::unit_x()); + let f = builder::tsweep(&e, Vector3::unit_y()); + let cube = builder::tsweep(&f, Vector3::unit_z()); + + let v = builder::vertex(Point3::new(0.5, 0.25, -0.5)); + let w = builder::rsweep(&v, Point3::new(0.5, 0.5, 0.0), Vector3::unit_z(), Rad(7.0)); + let f = builder::try_attach_plane(&[w]).unwrap(); + let mut cylinder = builder::tsweep(&f, Vector3::unit_z() * 2.0); + cylinder.not(); + let and = truck_shapeops::and(&cube, &cylinder, 0.05).unwrap(); + and.edge_iter().for_each(|edge| { + let mut curve = edge.curve(); + if let Curve::IntersectionCurve(inter) = &curve { + if matches! { inter.leader(), Leader::Polyline(_) } { + let flag = curve.to_bspline_leader(0.01, 0.1, 20); + println!("{flag}"); + } + } + edge.set_curve(curve); + }); + + let json = serde_json::to_vec_pretty(&and).unwrap(); + std::fs::write("punched-cube-shapeops.json", json).unwrap(); +} diff --git a/vendor/truck-shapeops/src/alternative.rs b/vendor/truck-shapeops/src/alternative.rs new file mode 100644 index 00000000..24af0b01 --- /dev/null +++ b/vendor/truck-shapeops/src/alternative.rs @@ -0,0 +1,182 @@ +use truck_meshalgo::prelude::*; + +#[derive(Clone, Debug)] +pub enum Alternative { + FirstType(T), + SecondType(U), +} + +#[doc(hidden)] +#[macro_export] +macro_rules! impl_from { + ($firsttype: ty, $secondtype: ty) => { + impl From<$firsttype> for $crate::alternative::Alternative<$firsttype, $secondtype> { + #[inline(always)] + fn from(t: $firsttype) -> Self { $crate::alternative::Alternative::FirstType(t) } + } + }; +} + +impl From for Alternative { + #[inline(always)] + fn from(u: U) -> Self { Alternative::SecondType(u) } +} + +// test for impl_from +impl_from!((), usize); + +macro_rules! derive_method { + ($method: tt, $return_type: ty, $($var: ident : $paramtype: ty),*) => { + fn $method (&self, $($var: $paramtype),*) -> $return_type { + match &self { + Alternative::FirstType(got) => got.$method($($var),*), + Alternative::SecondType(got) => got.$method($($var),*), + } + } + }; + ($method: tt <$x: ident : $y: path>, $return_type: ty, $($var: ident : $paramtype: ty),*) => { + fn $method<$x : $y>(&self, $($var: $paramtype),*) -> $return_type { + match &self { + Alternative::FirstType(got) => got.$method($($var),*), + Alternative::SecondType(got) => got.$method($($var),*), + } + } + }; +} + +impl ParametricCurve for Alternative +where + C0: ParametricCurve, + C1: ParametricCurve, +{ + type Point = C0::Point; + type Vector = C0::Vector; + derive_method!(subs, C0::Point, t: f64); + derive_method!(der, C0::Vector, t: f64); + derive_method!(der2, C0::Vector, t: f64); + derive_method!( + parameter_range, + (std::ops::Bound, std::ops::Bound), + ); +} + +impl BoundedCurve for Alternative +where + C0: BoundedCurve, + C1: BoundedCurve, +{ +} + +impl ParametricSurface for Alternative +where + S0: ParametricSurface, + S1: ParametricSurface, +{ + type Point = S0::Point; + type Vector = S0::Vector; + derive_method!(subs, S0::Point, u: f64, v: f64); + derive_method!(uder, S0::Vector, u: f64, v: f64); + derive_method!(vder, S0::Vector, u: f64, v: f64); + derive_method!(uuder, S0::Vector, u: f64, v: f64); + derive_method!(uvder, S0::Vector, u: f64, v: f64); + derive_method!(vvder, S0::Vector, u: f64, v: f64); +} + +impl ParametricSurface3D for Alternative +where + S0: ParametricSurface3D, + S1: ParametricSurface3D, +{ + derive_method!(normal, Vector3, u: f64, v: f64); +} + +impl Cut for Alternative +where + C0: Cut, + C1: Cut, +{ + fn cut(&mut self, t: f64) -> Self { + match self { + Self::FirstType(curve) => Self::FirstType(curve.cut(t)), + Self::SecondType(curve) => Self::SecondType(curve.cut(t)), + } + } +} + +impl ParameterDivision1D for Alternative +where + C0: ParameterDivision1D, + C1: ParameterDivision1D, +{ + type Point = C0::Point; + derive_method!( + parameter_division, + (Vec, Vec), + range: (f64, f64), + tol: f64 + ); +} + +impl ParameterDivision2D for Alternative +where + S0: ParameterDivision2D, + S1: ParameterDivision2D, +{ + derive_method!( + parameter_division, + (Vec, Vec), + range: ((f64, f64), (f64, f64)), + tol: f64 + ); +} + +impl SearchParameter for Alternative +where + T: SearchParameter, + U: SearchParameter, +{ + type Point = T::Point; + derive_method!( + search_parameter>, + Option, + point: T::Point, + hint: H, + trials: usize + ); +} + +impl SearchNearestParameter for Alternative +where + T: SearchNearestParameter, + U: SearchNearestParameter, +{ + type Point = T::Point; + derive_method!( + search_nearest_parameter>, + Option, + point: T::Point, + hint: H, + trials: usize + ); +} + +impl Invertible for Alternative +where + T: Invertible, + U: Invertible, +{ + #[inline(always)] + fn invert(&mut self) { + match self { + Self::FirstType(entity) => entity.invert(), + Self::SecondType(entity) => entity.invert(), + } + } + #[inline(always)] + fn inverse(&self) -> Self { + match self { + Self::FirstType(entity) => Self::FirstType(entity.inverse()), + Self::SecondType(entity) => Self::SecondType(entity.inverse()), + } + } +} diff --git a/vendor/truck-shapeops/src/healing/mod.rs b/vendor/truck-shapeops/src/healing/mod.rs new file mode 100644 index 00000000..df7b5e61 --- /dev/null +++ b/vendor/truck-shapeops/src/healing/mod.rs @@ -0,0 +1,132 @@ +use rustc_hash::{FxHashMap as HashMap, FxHashSet as HashSet}; +use truck_geometry::prelude::*; +use truck_meshalgo::rexport_polymesh::*; +use truck_topology::compress::*; + +type Edge = CompressedEdge; +type EdgeIndex = CompressedEdgeIndex; +type Wire = Vec; +type Face = CompressedFace; +type Shell = CompressedShell; + +trait SP: Fn(&S, Point3, Option<(f64, f64)>) -> Option<(f64, f64)> {} +impl SP for F where F: Fn(&S, Point3, Option<(f64, f64)>) -> Option<(f64, f64)> {} + +mod split_closed_edges; +use split_closed_edges::split_closed_edges; + +mod split_closed_faces; +use split_closed_faces::split_closed_faces; + +/// Splits closed edges and faces +/// +/// # Details +/// The topology of the shapes handled by truck has the following rules +/// - The endpoints of the edges must be different. +/// - The boundaries of the faces must be a simple wire. +/// +/// Shapes created in other CAD systems do not necessarily follow these rules. +/// When such shapes are handled by truck, this method is applied at the stage +/// of `CompressedShell` and `CompressedSolid`, which are intermediate forms. +/// +/// # Remarks +/// Boundary simplification is still only implemented for cylinders. +/// It has not yet been implemented for cases involving singularities, such as spherical surfaces. +pub trait SplitClosedEdgesAndFaces { + /// Splits closed edges and faces + fn split_closed_edges_and_faces(&mut self, tol: f64); +} + +impl SplitClosedEdgesAndFaces for CompressedShell +where + C: ParametricCurve3D + + BoundedCurve + + Cut + + ParameterDivision1D + + SearchNearestParameter + + TryFrom, S>>, + S: ParametricSurface3D + SearchParameter, +{ + fn split_closed_edges_and_faces(&mut self, tol: f64) { + fn sp(surface: &S, point: Point3, hint: Option<(f64, f64)>) -> Option<(f64, f64)> + where S: SearchParameter { + surface.search_parameter(point, hint, 100) + } + split_closed_edges(self); + split_closed_faces(self, tol, sp); + } +} + +impl SplitClosedEdgesAndFaces for CompressedSolid +where + C: ParametricCurve3D + + BoundedCurve + + Cut + + ParameterDivision1D + + SearchNearestParameter + + TryFrom, S>>, + S: ParametricSurface3D + SearchParameter, +{ + fn split_closed_edges_and_faces(&mut self, tol: f64) { + self.boundaries + .iter_mut() + .for_each(|shell| shell.split_closed_edges_and_faces(tol)) + } +} + +/// robust version of splitting closed edges and faces. +/// +/// # Details +/// Robust version of [`SplitClosedEdgesAndFaces`] based on [`SearchNearestParameter`]. +pub trait RobustSplitClosedEdgesAndFaces { + /// Splits closed edges and faces + fn robust_split_closed_edges_and_faces(&mut self, tol: f64); +} + +impl RobustSplitClosedEdgesAndFaces for CompressedShell +where + C: ParametricCurve3D + + BoundedCurve + + Cut + + ParameterDivision1D + + SearchNearestParameter + + TryFrom, S>>, + S: ParametricSurface3D + + SearchParameter + + SearchNearestParameter, +{ + fn robust_split_closed_edges_and_faces(&mut self, tol: f64) { + fn sp(surface: &S, point: Point3, hint: Option<(f64, f64)>) -> Option<(f64, f64)> + where S: SearchParameter + SearchNearestParameter + { + surface + .search_parameter(point, hint, 100) + .or_else(|| surface.search_parameter(point, None, 100)) + .or_else(|| surface.search_nearest_parameter(point, hint, 100)) + .or_else(|| surface.search_nearest_parameter(point, None, 100)) + } + split_closed_edges(self); + split_closed_faces(self, tol, sp); + } +} + +impl RobustSplitClosedEdgesAndFaces for CompressedSolid +where + C: ParametricCurve3D + + BoundedCurve + + Cut + + ParameterDivision1D + + SearchNearestParameter + + TryFrom, S>>, + S: ParametricSurface3D + + SearchParameter + + SearchNearestParameter, +{ + fn robust_split_closed_edges_and_faces(&mut self, tol: f64) { + let fs = RobustSplitClosedEdgesAndFaces::robust_split_closed_edges_and_faces; + self.boundaries.iter_mut().for_each(|shell| fs(shell, tol)) + } +} + +#[cfg(test)] +mod tests; diff --git a/vendor/truck-shapeops/src/healing/split_closed_edges.rs b/vendor/truck-shapeops/src/healing/split_closed_edges.rs new file mode 100644 index 00000000..666cd973 --- /dev/null +++ b/vendor/truck-shapeops/src/healing/split_closed_edges.rs @@ -0,0 +1,68 @@ +use super::*; + +/// Cuts the edge at the middle point of the bounded curve `curve`. +fn cut_edge( + vertices: &mut Vec

, + Edge { + vertices: (_, v1), + curve, + }: &mut Edge, +) -> Edge +where + C: BoundedCurve + Cut, +{ + let (t0, t1) = curve.range_tuple(); + let t = (t0 + t1) / 2.0; + vertices.push(curve.subs(t)); + + let v2 = *v1; + *v1 = vertices.len() - 1; + let v1 = *v1; + let curve0 = curve.cut(t); + Edge { + vertices: (v1, v2), + curve: curve0, + } +} + +/// Adds cut edges to the vector `edges`. +fn add_edges(vertices: &mut Vec

, edges: &mut Vec>) -> HashMap +where C: BoundedCurve + Cut { + let len = edges.len(); + let sub_add_edges = move |i: usize| { + let (v0, v1) = edges[i].vertices; + if v0 == v1 { + let new_edge = cut_edge(vertices, &mut edges[i]); + edges.push(new_edge); + Some((i, edges.len() - 1)) + } else { + None + } + }; + (0..len).filter_map(sub_add_edges).collect() +} + +fn replace_edges(wire: &mut Wire, added: &HashMap) { + let insert_one_edge = |edge: EdgeIndex| { + let Some(new_edge_index) = added.get(&edge.index) else { + return vec![edge]; + }; + let new_edge = EdgeIndex { + index: *new_edge_index, + orientation: edge.orientation, + }; + match edge.orientation { + true => vec![edge, new_edge], + false => vec![new_edge, edge], + } + }; + let new_wire = wire.iter().copied().flat_map(insert_one_edge).collect(); + *wire = new_wire; +} + +pub(super) fn split_closed_edges(shell: &mut Shell) +where C: BoundedCurve + Cut { + let added = add_edges(&mut shell.vertices, &mut shell.edges); + let wire_iter = shell.faces.iter_mut().flat_map(|face| &mut face.boundaries); + wire_iter.for_each(|wire| replace_edges(wire, &added)); +} diff --git a/vendor/truck-shapeops/src/healing/split_closed_faces.rs b/vendor/truck-shapeops/src/healing/split_closed_faces.rs new file mode 100644 index 00000000..41138830 --- /dev/null +++ b/vendor/truck-shapeops/src/healing/split_closed_faces.rs @@ -0,0 +1,975 @@ +use super::*; +use std::ops::Range; + +pub(super) fn split_closed_faces(shell: &mut Shell, tol: f64, sp: impl SP) +where + C: ParametricCurve3D + + BoundedCurve + + Cut + + ParameterDivision1D + + SearchNearestParameter + + TryFrom, S>>, + S: ParametricSurface3D, { + let to_poly = closure_to_poly(tol); + let mut poly_edges: Vec<_> = shell.edges.iter().map(to_poly).collect(); + let len = shell.faces.len(); + (0..len).for_each(|i| { + split_face_with_non_closed_boundary(i, shell, &mut poly_edges, &sp, tol); + }); + let len = shell.faces.len(); + let closure = |i| split_face_with_non_simple_wire(i, shell, &mut poly_edges, &sp, tol); + let new_faces: Vec<_> = (0..len).filter_map(closure).flatten().collect(); + shell.faces.extend(new_faces); +} + +fn split_face_with_non_closed_boundary( + face_index: usize, + shell: &mut Shell, + poly_edges: &mut Vec>, + sp: impl SP, + tol: f64, +) -> Option<()> +where + C: ParametricCurve3D + + BoundedCurve + + Cut + + ParameterDivision1D + + SearchNearestParameter + + TryFrom, S>>, + S: ParametricSurface3D, +{ + let (divisor, open, closed) = + non_closed_wires_in_param_divisor(face_index, shell, poly_edges, &sp)?; + debug_assert_eq!(open.len(), 2); + take_vertices_to_intersections(divisor, face_index, shell, poly_edges, &sp, tol); + let Shell { faces, edges, .. } = shell; + let new_boundaries = split_boundaries_by_divisor( + &faces[face_index], + &closed, + divisor, + edges, + poly_edges, + &sp, + tol, + )?; + let boundaries = &mut faces[face_index].boundaries; + connect_open_boundaries(boundaries, new_boundaries, divisor, edges, open) +} + +fn split_face_with_non_simple_wire( + face_index: usize, + shell: &mut Shell, + poly_edges: &mut Vec>, + sp: impl SP, + tol: f64, +) -> Option>> +where + C: ParametricCurve3D + + BoundedCurve + + Cut + + ParameterDivision1D + + SearchNearestParameter + + TryFrom, S>>, + S: ParametricSurface3D, +{ + let Face { boundaries, .. } = &shell.faces[face_index]; + let divisor = find_non_simple_wire_divisor(boundaries, &shell.edges, &shell.vertices)?; + split_face_by_divisor(face_index, divisor, shell, poly_edges, sp, tol) +} + +// --- find_non_simple_wire_divisor --- + +fn find_non_simple_wire_divisor( + boundaries: &[Wire], + edges: &[Edge], + vertices: &[Point3], +) -> Option<(usize, usize)> { + let closure = |boundary| non_simple_wire_divisor(boundary, edges, vertices); + boundaries.iter().cloned().find_map(closure) +} + +fn non_simple_wire_divisor( + mut boundary: Wire, + edges: &[Edge], + vertices: &[Point3], +) -> Option<(usize, usize)> { + let (i0, j0) = boundary.iter().enumerate().find_map(find_loop())?; + boundary.rotate_left(i0); + let j0 = j0 - i0; + + let (i1, j1) = boundary + .iter() + .enumerate() + .skip(j0 + 1) + .find_map(find_loop()) + .unwrap_or((j0, boundary.len())); + + let (k0, k1) = ((j0 + 1) / 2, (i1 + j1 + 1) / 2); + let f = closure_take_front(edges); + let pre_divisor = (f(boundary[k0]), f(boundary[k1])); + nearest_correction( + pre_divisor, + &boundary[1..j0], + &boundary[i1 + 1..j1], + vertices, + closure_take_front(edges), + ) +} + +fn find_loop() -> impl FnMut((usize, &EdgeIndex)) -> Option<(usize, usize)> { + let mut map = HashMap::::default(); + move |(i, edge)| map.insert(edge.index, i).map(move |j| (j, i)) +} + +// --- find_non_closed_wires_in_param_divisor --- + +type FindNonClosedWiresInParamDivisorResult = Option<((usize, usize), Vec, Vec)>; +fn non_closed_wires_in_param_divisor( + face_index: usize, + Shell { + edges, + faces, + ref vertices, + }: &mut Shell, + poly_edges: &[PolylineCurve], + sp: impl SP, +) -> FindNonClosedWiresInParamDivisorResult +where + C: ParametricCurve3D + BoundedCurve + TryFrom, S>>, + S: ParametricSurface3D, +{ + let Face { + boundaries, + surface, + .. + } = &faces[face_index]; + let param_boundaries = create_param_boundaries(boundaries, surface, poly_edges, sp)?; + let (mut open, mut closed) = (Vec::new(), Vec::new()); + boundaries + .iter() + .zip(¶m_boundaries) + .for_each(|(boundary, param_boundary)| { + let poly = PolylineCurve(param_boundary.iter().flatten().copied().collect()); + match boundary_type(&poly) { + BoundaryType::NotClosed => open.push((boundary, param_boundary)), + _ => closed.push((boundary, param_boundary)), + } + }); + if open.len() != 2 { + return None; + } + let take_front = &closure_take_front(edges); + let find_nearest_to_ends = move |i: usize, u0: f64, up: f64| { + move |(edge_index, poly): ZippedEdge<'_>| { + let u = poly[0][i]; + let closure = move |i| { + let u = u + i as f64 * up; + f64::abs(u - u0) + }; + let dist = (-2..=2) + .map(closure) + .min_by(|x, y| x.partial_cmp(y).unwrap()) + .unwrap(); + (take_front(*edge_index), dist) + } + }; + let create_pre_divisor = move |(wire, poly_wire): ZippedWire<'_>| { + let (urange, vrange) = surface.try_range_tuple(); + let (up, vp) = (surface.u_period(), surface.v_period()); + let p = poly_wire[0][0]; + let q = *poly_wire.last().unwrap().last().unwrap(); + let vertices: Vec<(usize, f64)> = if p.x.near(&q.x) { + if let (Some(vp), Some((v0, _))) = (vp, vrange) { + let closure = find_nearest_to_ends(1, v0, vp); + Some(wire.iter().zip(poly_wire).map(closure).collect()) + } else { + None + } + } else if p.y.near(&q.y) { + if let (Some(up), Some((u0, _))) = (up, urange) { + let closure = find_nearest_to_ends(0, u0, up); + Some(wire.iter().zip(poly_wire).map(closure).collect()) + } else { + None + } + } else { + None + }?; + vertices + .into_iter() + .min_by(|(_, x), (_, y)| x.partial_cmp(y).unwrap()) + .map(|(v, _)| v) + }; + let pre_divisor_vec = open + .iter() + .copied() + .map(create_pre_divisor) + .collect::>>()?; + let open: Vec<_> = open.into_iter().map(|(w, _)| w.to_vec()).collect(); + let closed: Vec<_> = closed.into_iter().map(|(w, _)| w.to_vec()).collect(); + let divisor = nearest_correction( + (pre_divisor_vec[0], pre_divisor_vec[1]), + &open[0], + &open[1], + vertices, + take_front, + )?; + Some((divisor, open, closed)) +} + +// --- split_face_by_divisor --- + +fn split_face_by_divisor( + face_index: usize, + divisor: (usize, usize), + shell: &mut Shell, + poly_edges: &mut Vec>, + sp: impl SP, + tol: f64, +) -> Option>> +where + C: ParametricCurve3D + + BoundedCurve + + Cut + + ParameterDivision1D + + SearchNearestParameter + + TryFrom, S>>, + S: ParametricSurface3D, +{ + take_vertices_to_intersections(divisor, face_index, shell, poly_edges, &sp, tol); + let face = &shell.faces[face_index]; + let new_boundaries = split_boundaries_by_divisor( + face, + &face.boundaries, + divisor, + &mut shell.edges, + poly_edges, + &sp, + tol, + )?; + divide_face( + &mut shell.faces[face_index], + new_boundaries, + poly_edges, + &sp, + ) +} + +// --- take_vertices_to_intersections --- + +fn take_vertices_to_intersections( + divisor: (usize, usize), + face_index: usize, + Shell { + vertices, + edges, + faces, + }: &mut Shell, + poly_edges: &mut Vec>, + sp: impl SP, + tol: f64, +) -> Option<()> +where + C: ParametricCurve3D + + BoundedCurve + + Cut + + SearchNearestParameter + + ParameterDivision1D, + S: ParametricSurface3D, +{ + let Face { + ref boundaries, + ref surface, + .. + } = faces[face_index]; + let param_boundaries = create_param_boundaries(boundaries, surface, poly_edges, &sp)?; + let (v0, v1) = divisor; + let p = zip_boundaries(boundaries, ¶m_boundaries) + .find_map(closure_find_vertex_parameter(v0, edges))?; + let q = zip_boundaries(boundaries, ¶m_boundaries) + .find_map(closure_find_vertex_parameter(v1, edges))?; + let Face { + ref boundaries, + ref surface, + .. + } = faces[face_index]; + let periods = (surface.u_period(), surface.v_period()); + let q = periodic_iterator(q, periods) + .min_by(|q, r| p.distance2(*q).partial_cmp(&p.distance2(*r)).unwrap())?; + let pcurve = PCurve::new(Line(p, q), surface.clone()); + let cut_edge = |(edge_index, param_edge): ZippedEdge<'_>| { + let index = edge_index.index; + let vec = enumerate_intersections(&edges[index], param_edge, &pcurve)?; + if vec.is_empty() { + return Some(None); + } + let erange = cut_edge_by_intersections(index, vec, edges, poly_edges, vertices, tol); + Some(Some((index, erange))) + }; + let new_edges = zip_boundaries(boundaries, ¶m_boundaries) + .map(cut_edge) + .collect::>>()?; + let insert = |(index, erange): (usize, Range)| { + faces + .iter_mut() + .flat_map(|face| &mut face.boundaries) + .for_each(|wire| insert_new_edges(wire, index, erange.clone())); + }; + new_edges.into_iter().flatten().for_each(insert); + Some(()) +} + +fn closure_find_vertex_parameter( + v: usize, + edges: &[Edge], +) -> impl Fn(ZippedEdge<'_>) -> Option + '_ { + let take_front = closure_take_front(edges); + move |(edge_index, param_edge): ZippedEdge<'_>| match v == take_front(*edge_index) { + true => Some(param_edge[0]), + false => None, + } +} + +fn enumerate_intersections( + edge: &Edge, + param_edge: &PolylineCurve, + pcurve: &PCurve, S>, +) -> Option> +where + C: ParametricCurve3D + SearchNearestParameter, + S: ParametricSurface3D, +{ + let rough_intersections = intersections_between_line_polyline(*pcurve.curve(), param_edge); + let intersections = exact_intersections(rough_intersections, pcurve, &edge.curve)?; + Some(intersections) +} + +fn intersections_between_line_polyline( + line: Line, + param_edge: &PolylineCurve, +) -> Vec { + let filter = |p: &[Point2]| { + let (s, t, _) = line.intersection(Line(p[0], p[1]))?; + let unit = 0.0..1.0; + match unit.contains(&s) && unit.contains(&t) { + true => Some(s), + false => None, + } + }; + let mut vec: Vec<_> = param_edge.windows(2).filter_map(filter).collect(); + if !vec.is_empty() { + let first = *param_edge.first().unwrap(); + let last = *param_edge.last().unwrap(); + if line.distance_to_point(first).so_small() { + vec.remove(0); + } + if line.distance_to_point(last).so_small() { + vec.pop(); + } + } + vec +} + +fn exact_intersections( + naive_intersections: Vec, + pcurve: &PCurve, S>, + curve: &C, +) -> Option> +where + C: ParametricCurve3D + SearchNearestParameter, + S: ParametricSurface3D, +{ + let closure = |t| { + let (_, t, p) = search_intersection(pcurve, curve, t)?; + Some((t, p)) + }; + let mut res = naive_intersections + .into_iter() + .map(closure) + .collect::>>()?; + res.sort_by(|(s, _), (t, _)| s.partial_cmp(t).unwrap()); + Some(res) +} + +fn search_intersection( + pcurve: &PCurve, S>, + curve: &C, + mut t0: f64, +) -> Option<(f64, f64, Point3)> +where + C: ParametricCurve3D + SearchNearestParameter, + S: ParametricSurface3D, +{ + let mut previous0 = t0; + let mut previous1 = None; + for _ in 0..100 { + let p = pcurve.subs(t0); + let t1 = curve.search_nearest_parameter(p, previous1, 100)?; + let q = curve.subs(t1); + t0 = pcurve.search_nearest_parameter(q, t0, 100)?; + if let Some(previous1) = previous1 { + if previous0.near(&t0) && previous1.near(&t1) { + return Some((t0, t1, pcurve.subs(t0))); + } + } + previous0 = t0; + previous1 = Some(t1); + } + None +} + +fn cut_edge_by_intersections( + edge_index: usize, + intersections: Vec<(f64, Point3)>, + edges: &mut Vec>, + poly_edges: &mut Vec>, + vertices: &mut Vec, + tol: f64, +) -> Range +where + C: Cut + ParameterDivision1D, +{ + let edge = &mut edges[edge_index]; + let vfirst = vertices.len(); + vertices.extend(intersections.iter().map(|(_, p)| *p)); + let intersections_iter = intersections.into_iter().enumerate().rev(); + let closure = |(i, (t, _))| Edge { + vertices: (vfirst + i, vfirst + i + 1), + curve: edge.curve.cut(t), + }; + let mut new_edges: Vec<_> = intersections_iter.map(closure).collect(); + let to_poly = closure_to_poly(tol); + poly_edges[edge_index] = to_poly(edge); + poly_edges.extend(new_edges.iter().map(to_poly)); + new_edges.last_mut().unwrap().vertices.1 = edge.vertices.1; + edge.vertices.1 = vfirst; + let efirst = edges.len(); + edges.extend(new_edges); + efirst..edges.len() +} + +fn insert_new_edges(wire: &mut Wire, pivot_edge_index: usize, inserted_range: Range) { + let positions = wire + .iter() + .enumerate() + .filter_map( + |(i, edge_index)| match edge_index.index == pivot_edge_index { + true => Some(i), + false => None, + }, + ) + .collect::>(); + positions.into_iter().rev().for_each(|i| { + let orientation = wire[i].orientation; + let ei = |index| EdgeIndex { index, orientation }; + if orientation { + let tmp_wire = wire.split_off(i + 1); + wire.extend(inserted_range.clone().map(ei)); + wire.extend(tmp_wire); + } else { + let tmp_wire = wire.split_off(i); + wire.extend(inserted_range.clone().rev().map(ei)); + wire.extend(tmp_wire); + } + }); +} + +// --- split boundaries --- + +fn split_boundaries_by_divisor( + Face { + boundaries, + surface, + .. + }: &Face, + closed: &[Wire], + divisor: (usize, usize), + edges: &mut Vec>, + poly_edges: &mut Vec>, + sp: impl SP, + tol: f64, +) -> Option> +where + C: ParametricCurve3D + BoundedCurve + TryFrom, S>>, + S: ParametricSurface3D, +{ + let param_boundaries = create_param_boundaries(boundaries, surface, poly_edges, &sp)?; + let param_vertices = create_param_vertices(boundaries, ¶m_boundaries, edges); + let mut duplicated_edges = duplicated_edges(boundaries.iter().flatten().map(|ei| ei.index)); + let vertices_on_divisor = enumerate_vertices_on_divisor(divisor, ¶m_vertices, surface)?; + let new_edges = create_new_edges(&vertices_on_divisor, poly_edges, surface, tol)?; + let (new_edge_range, new_edge_indices) = signup_new_edges(edges, new_edges); + duplicated_edges.extend(new_edge_range); + let edge_iter = closed.iter().flatten().copied().chain(new_edge_indices); + let vemap = create_vemap(edges, edge_iter); + Some(construct_boundaries(vemap, edges, &duplicated_edges)) +} + +// --- create_param_boundaries --- + +fn create_param_boundaries( + boundaries: &[Wire], + surface: &S, + poly_edges: &[PolylineCurve], + sp: impl SP, +) -> Option>>> { + let closure = |boundary: &Wire| boundary_to_param_polys(boundary, surface, poly_edges, &sp); + boundaries.iter().map(closure).collect() +} + +fn boundary_to_param_polys( + boundary: &Wire, + surface: &S, + poly_edges: &[PolylineCurve], + sp: impl SP, +) -> Option>> { + let get_poly = closure_get_poly(poly_edges); + let poly_boundary: Vec<_> = boundary.iter().copied().map(get_poly).collect(); + let poly_wire_iter = PolyWireIter::try_new(&poly_boundary)?; + let closure = poly_project_to_uv(surface, sp); + let mut long_poly = poly_wire_iter.map(closure).collect::>>()?; + boundary_into_domain(&mut long_poly, surface); + let mut vec = vec![long_poly.remove(0)]; + let split_into_2dpoly = move |poly: &PolylineCurve<_>| { + let latter = long_poly.split_off(poly.len() - 1); + vec.append(&mut long_poly); + let mut res = Vec::new(); + res.append(&mut vec); + vec.push(*res.last().unwrap()); + long_poly = latter; + PolylineCurve(res) + }; + Some(poly_boundary.iter().map(split_into_2dpoly).collect()) +} + +struct PolyWireIter<'a, P> { + current: std::slice::Iter<'a, P>, + stock: std::vec::IntoIter>, +} + +impl<'a, P> PolyWireIter<'a, P> { + fn try_new + 'a>(vec: &'a [I]) -> Option { + let iters: Vec<_> = vec.iter().map(|v| v.as_ref().iter()).collect(); + let mut stock = iters.into_iter(); + Some(Self { + current: stock.next()?, + stock, + }) + } +} + +impl<'a, P: Copy> Iterator for PolyWireIter<'a, P> { + type Item = P; + fn next(&mut self) -> Option { + if let Some(next) = self.current.next() { + Some(*next) + } else { + self.current = self.stock.next()?; + self.current.next()?; + self.current.next().cloned() + } + } +} + +fn poly_project_to_uv<'a, S: ParametricSurface3D>( + surface: &'a S, + sp: impl SP + 'a, +) -> impl FnMut(Point3) -> Option + 'a { + let (up, vp) = (surface.u_period(), surface.v_period()); + let mut previous = None; + move |pt| { + let (mut u, mut v) = sp(surface, pt, previous)?; + if let (Some(up), Some((u0, _))) = (up, previous) { + u = get_mindiff(u, u0, up); + } + if let (Some(vp), Some((_, v0))) = (vp, previous) { + v = get_mindiff(v, v0, vp); + } + previous = Some((u, v)); + Some(Point2::new(u, v)) + } +} + +fn abs_diff(previous: f64) -> impl Fn(&f64, &f64) -> std::cmp::Ordering { + let f = move |x: &f64| f64::abs(x - previous); + move |x: &f64, y: &f64| f(x).partial_cmp(&f(y)).unwrap() +} +fn get_mindiff(u: f64, u0: f64, up: f64) -> f64 { + let closure = |i| u + i as f64 * up; + (-2..=2).map(closure).min_by(abs_diff(u0)).unwrap() +} + +fn boundary_into_domain(vec: &mut Vec, surface: &S) { + let (up, vp) = (surface.u_period(), surface.v_period()); + let (urange, vrange) = surface.try_range_tuple(); + let grav = vec.iter().fold(Point2::origin(), |g, p| g + p.to_vec()) / vec.len() as f64; + if let (Some(up), Some((u0, _))) = (up, urange) { + let quot = f64::floor((grav.x - u0) / up); + vec.iter_mut().for_each(|p| p.x -= quot * up); + } + if let (Some(vp), Some((v0, _))) = (vp, vrange) { + let quot = f64::floor((grav.y - v0) / vp); + vec.iter_mut().for_each(|p| p.y -= quot * vp); + } + let last = *vec.last().unwrap(); + if !vec[0].near(&last) { + let Point2 { x: u0, y: v0 } = last; + if surface.uder(u0, v0).so_small() || surface.vder(u0, v0).so_small() { + vec.push(vec[0]); + } + } +} + +// --- create_param_vertices --- + +fn create_param_vertices( + boundaries: &[Wire], + param_boundaries: &[Vec>], + edges: &[Edge], +) -> HashMap { + let take_front = closure_take_front(edges); + zip_boundaries(boundaries, param_boundaries) + .map(|(edge_index, param_edge)| (take_front(*edge_index), param_edge[0])) + .collect() +} + +// --- duplicated_edges --- + +fn duplicated_edges(edges: impl Iterator) -> HashSet { + let mut done = HashSet::default(); + edges.filter(move |index| !done.insert(*index)).collect() +} + +// --- enumerate_vertices_on_divisor --- + +fn enumerate_vertices_on_divisor( + divisor: (usize, usize), + param_vertices: &HashMap, + surface: &S, +) -> Option> { + let (v0, v1) = divisor; + let (p, q) = (*param_vertices.get(&v0)?, *param_vertices.get(&v1)?); + let periods = (surface.u_period(), surface.v_period()); + let q = periodic_iterator(q, periods) + .min_by(|q, r| p.distance2(*q).partial_cmp(&p.distance2(*r)).unwrap())?; + let line = Line(p, q); + let iter = param_vertices.iter().filter_map(move |(v, uv)| { + periodic_iterator(*uv, periods) + .find(move |uv| line.distance_to_point_as_segment(*uv).so_small()) + .map(move |uv| Some((line.search_nearest_parameter(uv, None, 1)?, (*v, uv)))) + }); + let mut vertices_on_divisor = iter.collect::>>()?; + vertices_on_divisor.sort_by(|(s, _), (t, _)| s.partial_cmp(t).unwrap()); + Some(vertices_on_divisor) +} + +// --- create_new_edges --- + +fn create_new_edges( + vertices_on_divisor: &[(f64, (usize, Point2))], + poly_edges: &mut Vec>, + surface: &S, + tol: f64, +) -> Option>> +where + C: TryFrom, S>>, + S: ParametricSurface3D, +{ + vertices_on_divisor + .iter() + .for_each(|(_, (_, uv))| print!("{:?} ", surface.subs(uv.x, uv.y))); + println!(); + let make_edge = move |p: &[(f64, (usize, Point2))]| { + let ((_, (v0, uv0)), (_, (v1, uv1))) = (p[0], p[1]); + let pcurve = PCurve::new(Line(uv0, uv1), surface.clone()); + poly_edges.push(PolylineCurve::from_curve(&pcurve, (0.0, 1.0), tol)); + Some(Edge { + vertices: (v0, v1), + curve: C::try_from(pcurve).ok()?, + }) + }; + vertices_on_divisor.chunks(2).map(make_edge).collect() +} + +// --- signup new edges --- + +fn signup_new_edges( + edges: &mut Vec>, + new_edges: Vec>, +) -> (Range, impl Iterator) { + let len = edges.len(); + edges.extend(new_edges); + let new_edge_indices = (len..edges.len()).flat_map(move |index| { + let ei = |orientation: bool| EdgeIndex { index, orientation }; + [ei(true), ei(false)] + }); + (len..edges.len(), new_edge_indices) +} + +// --- create_vemap --- + +fn create_vemap( + edges: &[Edge], + edge_iter: impl Iterator, +) -> HashMap> { + let mut vemap = HashMap::>::default(); + let take_front = closure_take_front(edges); + edge_iter.for_each(|edge_index| { + let v = take_front(edge_index); + vemap.entry(v).or_insert_with(Vec::new).push(edge_index); + }); + vemap +} + +// --- construct boundaries --- + +fn construct_boundaries( + mut vemap: HashMap>, + edges: &[Edge], + new_edge_range: &HashSet, +) -> Vec { + let take_back = closure_take_back(edges); + let mut new_boundaries = Vec::new(); + while !vemap.is_empty() { + let (start, vec) = vemap.iter_mut().min_by_key(|(idx, _)| *idx).unwrap(); + let start = *start; + let mut edge_index = vec.pop().unwrap(); + if vec.is_empty() { + vemap.remove(&start); + } + let mut wire = Vec::new(); + loop { + wire.push(edge_index); + let v = take_back(edge_index); + if v == start { + break; + } + let Some(vec) = vemap.get_mut(&v) else { + unreachable!(); + }; + match vec.len() { + 1 => { + edge_index = vec.pop().unwrap(); + vemap.remove(&v); + } + 2 => { + let is_new0 = new_edge_range.contains(&edge_index.index); + let is_new1 = new_edge_range.contains(&vec[0].index); + let i = if is_new0 != is_new1 { 0 } else { 1 }; + edge_index = vec.remove(i); + } + _ => panic!("something wrong!"), + } + } + new_boundaries.push(wire); + } + new_boundaries +} + +// --- divide_face --- + +fn divide_face( + Face { + ref mut boundaries, + ref surface, + ref orientation, + }: &mut Face, + new_boundaries: Vec, + poly_edges: &[PolylineCurve], + sp: impl SP, +) -> Option>> { + let mut face_boundaries = assort_boundary(surface, new_boundaries, poly_edges, sp)?.into_iter(); + *boundaries = face_boundaries.next().unwrap(); + let create_face = |boundaries| Face { + boundaries, + surface: surface.clone(), + orientation: *orientation, + }; + Some(face_boundaries.map(create_face).collect()) +} + +fn assort_boundary( + surface: &S, + boundaries: Vec, + poly_edges: &[PolylineCurve], + sp: impl SP, +) -> Option>> +where + S: ParametricSurface3D, +{ + let (mut positives, mut negatives) = (Vec::new(), Vec::new()); + let param_boundaries = create_param_boundaries(&boundaries, surface, poly_edges, sp)?; + boundaries + .into_iter() + .zip(param_boundaries) + .try_for_each(|(boundary, param_boundary)| { + let poly_boundary = PolylineCurve(param_boundary.into_iter().flatten().collect()); + match boundary_type(&poly_boundary) { + BoundaryType::Positive => positives.push((boundary, poly_boundary)), + BoundaryType::Negative => negatives.push((boundary, poly_boundary)), + _ => return None, + } + Some(()) + })?; + let mut assorted = positives.into_iter().map(|x| vec![x]).collect::>(); + negatives + .into_iter() + .try_for_each(|(boundary, poly_boundary)| { + let p = poly_boundary[0]; + let vec = assorted.iter_mut().find(|vec| vec[0].1.include(p))?; + vec.push((boundary, poly_boundary)); + Some(()) + })?; + let res = assorted + .into_iter() + .map(|vec| vec.into_iter().map(|(x, _)| x).collect()) + .collect(); + Some(res) +} + +// --- connect_open_boundaries --- + +fn connect_open_boundaries( + boundaries: &mut Vec, + mut new_boundaries: Vec, + divisor: (usize, usize), + edges: &[Edge], + mut open: Vec, +) -> Option<()> { + let (v0, v1) = divisor; + let take_front = closure_take_front(edges); + let find_index = move |wire: &Wire, v: usize| { + wire.iter() + .position(|edge_index| take_front(*edge_index) == v) + }; + let (i, index) = open + .iter() + .enumerate() + .find_map(|(i, wire)| find_index(wire, v0).map(|index| (i, index)))?; + if index != 0 { + open.swap(0, 1); + } + let (mut open_v1, mut open_v0) = (open.pop()?, open.pop()?); + let j = find_index(&open_v1, v1)?; + open_v0.rotate_left(i); + open_v1.rotate_left(j); + let (k, added_wire) = new_boundaries + .iter_mut() + .find_map(|wire| find_index(wire, v0).map(|index| (index, wire)))?; + added_wire.rotate_left(k); + open_v0.append(added_wire); + let l = find_index(&open_v0, v1)?; + open_v0.rotate_left(l); + open_v0.append(&mut open_v1); + *added_wire = open_v0; + *boundaries = new_boundaries; + Some(()) +} + +// --- common functions --- + +fn closure_take_front(edges: &[Edge]) -> impl Fn(EdgeIndex) -> usize + '_ { + move |edge_index: EdgeIndex| { + let index = edge_index.index; + let (v0, v1) = edges[index].vertices; + match edge_index.orientation { + true => v0, + false => v1, + } + } +} + +fn closure_take_back(edges: &[Edge]) -> impl Fn(EdgeIndex) -> usize + '_ { + move |edge_index: EdgeIndex| { + let (v0, v1) = edges[edge_index.index].vertices; + match edge_index.orientation { + true => v1, + false => v0, + } + } +} + +fn closure_get_poly( + poly_edges: &[PolylineCurve

], +) -> impl Fn(EdgeIndex) -> PolylineCurve

+ '_ { + move |edge_index: EdgeIndex| { + let poly = &poly_edges[edge_index.index]; + match edge_index.orientation { + true => poly.clone(), + false => poly.inverse(), + } + } +} + +fn closure_to_poly(tol: f64) -> impl Fn(&Edge) -> PolylineCurve

+ 'static +where C: BoundedCurve + ParameterDivision1D { + move |Edge { curve, .. }: &Edge| PolylineCurve::from_curve(curve, curve.range_tuple(), tol) +} + +type ZippedEdge<'a> = (&'a EdgeIndex, &'a PolylineCurve); +type ZippedWire<'a> = (&'a Wire, &'a Vec>); + +fn zip_boundaries<'a>( + boundaries: &'a [Wire], + param_boundaries: &'a [Vec>], +) -> impl Iterator> { + boundaries + .iter() + .zip(param_boundaries) + .flat_map(move |(b, pb)| b.iter().zip(pb)) +} + +#[derive(Clone, Copy, Debug)] +enum BoundaryType { + Positive, + Negative, + NotClosed, +} + +fn boundary_type(poly: &PolylineCurve) -> BoundaryType { + if poly[0].distance2(poly[poly.len() - 1]) < 1.0e-3 { + match poly.area() > 0.0 { + true => BoundaryType::Positive, + false => BoundaryType::Negative, + } + } else { + BoundaryType::NotClosed + } +} + +fn periodic_iterator( + p: Point2, + (up, vp): (Option, Option), +) -> impl Iterator { + let up_range = match up { + Some(up) => vec![-2.0 * up, -up, 0.0, up, 2.0 * up], + None => vec![0.0], + }; + let vp_range = match vp { + Some(vp) => vec![-2.0 * vp, -vp, 0.0, vp, 2.0 * vp], + None => vec![0.0], + }; + itertools::iproduct!(up_range, vp_range).map(move |(dx, dy)| p + Vector2::new(dx, dy)) +} + +fn nearest_correction( + (v0, _v1): (usize, usize), + wire0: &[EdgeIndex], + wire1: &[EdgeIndex], + vertices: &[Point3], + take_front: impl Fn(EdgeIndex) -> usize, +) -> Option<(usize, usize)> { + let p0 = vertices[v0]; + let v1 = wire1.iter().copied().map(&take_front).min_by(|v, w| { + p0.distance2(vertices[*v]) + .partial_cmp(&p0.distance2(vertices[*w])) + .unwrap() + })?; + let p1 = vertices[v1]; + let v0 = wire0.iter().copied().map(&take_front).min_by(|v, w| { + p1.distance2(vertices[*v]) + .partial_cmp(&p1.distance2(vertices[*w])) + .unwrap() + })?; + Some((v0, v1)) +} diff --git a/vendor/truck-shapeops/src/healing/tests.rs b/vendor/truck-shapeops/src/healing/tests.rs new file mode 100644 index 00000000..deb10588 --- /dev/null +++ b/vendor/truck-shapeops/src/healing/tests.rs @@ -0,0 +1,1534 @@ +use super::*; +use std::f64::consts::PI; +use truck_topology::Shell; + +#[test] +fn test_split_closed_edges() { + let vertices = vec![Point2::new(1.0, 0.0)]; + let curve = TrimmedCurve::new(UnitCircle::::new(), (0.0, 2.0 * PI)); + let edges = vec![CompressedEdge { + vertices: (0, 0), + curve, + }]; + let faces = vec![ + CompressedFace { + surface: (), + orientation: true, + boundaries: vec![vec![CompressedEdgeIndex { + index: 0, + orientation: true, + }]], + }, + CompressedFace { + surface: (), + orientation: false, + boundaries: vec![vec![CompressedEdgeIndex { + index: 0, + orientation: false, + }]], + }, + ]; + let mut shell = CompressedShell { + vertices, + edges, + faces, + }; + assert!(Shell::extract(shell.clone()).is_err()); + + split_closed_edges(&mut shell); + assert!(Shell::extract(shell.clone()).is_ok()); + + let CompressedShell { + vertices, + edges, + faces, + } = &shell; + assert_eq!(vertices.len(), 2); + assert_near!(vertices[0], Point2::new(1.0, 0.0)); + assert_near!(vertices[1], Point2::new(-1.0, 0.0)); + + assert_eq!(edges.len(), 2); + assert_eq!(edges[0].vertices, (0, 1)); + assert_eq!(edges[1].vertices, (1, 0)); + + assert_eq!( + *faces, + vec![ + CompressedFace { + surface: (), + orientation: true, + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 0, + orientation: true, + }, + CompressedEdgeIndex { + index: 1, + orientation: true, + } + ]], + }, + CompressedFace { + surface: (), + orientation: false, + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 1, + orientation: false, + }, + CompressedEdgeIndex { + index: 0, + orientation: false, + } + ]], + }, + ] + ); +} + +#[test] +fn test_split_closed_face_simple_cylinder_case() { + type Surface = RevolutedCurve>; + #[derive( + Clone, + Debug, + ParametricCurve, + BoundedCurve, + ParameterDivision1D, + Cut, + SearchNearestParameterD1, + )] + enum Curve { + Line(Line), + Arc(TrimmedCurve, Matrix4>>), + #[allow(clippy::enum_variant_names)] + PCurve(PCurve, Surface>), + } + impl From, Surface>> for Curve { + fn from(value: PCurve, Surface>) -> Self { Self::PCurve(value) } + } + + let vertices = vec![ + Point3::new(1.0, 0.0, 0.0), + Point3::new(-1.0, 0.0, 0.0), + Point3::new(1.0, 0.0, 1.0), + Point3::new(-1.0, 0.0, 1.0), + ]; + + let translate = Matrix4::from_translation(Vector3::unit_z()); + let edges = vec![ + CompressedEdge { + vertices: (0, 1), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()), + (0.0, PI), + )), + }, + CompressedEdge { + vertices: (1, 0), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()), + (PI, 2.0 * PI), + )), + }, + CompressedEdge { + vertices: (0, 2), + curve: Curve::Line(Line(vertices[0], vertices[2])), + }, + CompressedEdge { + vertices: (2, 3), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed(translate), + (0.0, PI), + )), + }, + CompressedEdge { + vertices: (3, 2), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed(translate), + (PI, 2.0 * PI), + )), + }, + ]; + let surface = RevolutedCurve::by_revolution( + Line(vertices[2], vertices[0]), + Point3::origin(), + Vector3::unit_z(), + ); + let faces = vec![Face { + surface, + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 1, + orientation: true, + }, + CompressedEdgeIndex { + index: 2, + orientation: true, + }, + CompressedEdgeIndex { + index: 4, + orientation: false, + }, + CompressedEdgeIndex { + index: 3, + orientation: false, + }, + CompressedEdgeIndex { + index: 2, + orientation: false, + }, + CompressedEdgeIndex { + index: 0, + orientation: true, + }, + ]], + orientation: true, + }]; + let mut shell = CompressedShell { + vertices, + edges, + faces, + }; + + assert!(Shell::extract(shell.clone()).is_err()); + + split_closed_faces(&mut shell, 0.01, sp); + + assert!(Shell::extract(shell.clone()).is_ok()); + let CompressedShell { + ref vertices, + ref edges, + ref mut faces, + } = shell; + assert_eq!(vertices.len(), 4); + assert_eq!(edges.len(), 6); + assert_eq!(edges[5].vertices, (3, 1)); + let curve0 = &edges[5].curve; + let curve1 = Line(Point3::new(-1.0, 0.0, 1.0), Point3::new(-1.0, 0.0, 0.0)); + for i in 0..=10 { + let t = i as f64 / 10.0; + assert_near!(curve0.subs(t), curve1.subs(t)); + } + assert_eq!(faces.len(), 2); + let i = faces + .iter_mut() + .position(|face| { + face.boundaries[0].contains(&CompressedEdgeIndex { + index: 2, + orientation: true, + }) + }) + .unwrap(); + if i == 1 { + faces.swap(0, 1); + } + let i = faces[0].boundaries[0] + .iter() + .position(|edge_index| { + *edge_index + == CompressedEdgeIndex { + index: 2, + orientation: true, + } + }) + .unwrap(); + faces[0].boundaries[0].rotate_left(i); + let i = faces[1].boundaries[0] + .iter() + .position(|edge_index| { + *edge_index + == CompressedEdgeIndex { + index: 3, + orientation: false, + } + }) + .unwrap(); + faces[1].boundaries[0].rotate_left(i); + + assert_eq!( + *faces, + vec![ + Face { + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 2, + orientation: true, + }, + CompressedEdgeIndex { + index: 4, + orientation: false, + }, + CompressedEdgeIndex { + index: 5, + orientation: true, + }, + CompressedEdgeIndex { + index: 1, + orientation: true, + } + ]], + surface, + orientation: true, + }, + Face { + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 3, + orientation: false, + }, + CompressedEdgeIndex { + index: 2, + orientation: false, + }, + CompressedEdgeIndex { + index: 0, + orientation: true, + }, + CompressedEdgeIndex { + index: 5, + orientation: false, + }, + ]], + surface, + orientation: true, + }, + ] + ); +} + +#[test] +fn test_split_closed_face_cylinder_with_hole() { + type Surface = RevolutedCurve>; + #[derive( + Clone, + Debug, + ParametricCurve, + BoundedCurve, + ParameterDivision1D, + Cut, + SearchNearestParameterD1, + )] + enum ParameterCurve { + Line(Line), + Arc(TrimmedCurve, Matrix3>>), + } + #[derive( + Clone, + Debug, + ParametricCurve, + BoundedCurve, + ParameterDivision1D, + Cut, + SearchNearestParameterD1, + )] + enum Curve { + Line(Line), + Arc(TrimmedCurve, Matrix4>>), + PCurve(PCurve), + } + impl From, Surface>> for Curve { + fn from(value: PCurve, Surface>) -> Self { + let (line, surface) = value.decompose(); + Self::PCurve(PCurve::new(ParameterCurve::Line(line), surface)) + } + } + + let vertices = vec![ + Point3::new(1.0, 0.0, 0.0), + Point3::new(-1.0, 0.0, 0.0), + Point3::new(1.0, 0.0, 1.0), + Point3::new(-1.0, 0.0, 1.0), + Point3::new(-1.0, 0.0, 0.25), + Point3::new(-1.0, 0.0, 0.75), + ]; + + let surface = RevolutedCurve::by_revolution( + Line(vertices[2], vertices[0]), + Point3::origin(), + Vector3::unit_z(), + ); + + let translate = Matrix4::from_translation(Vector3::unit_z()); + let transform = Matrix3::from_translation(Vector2::new(0.5, PI)) * Matrix3::from_scale(0.25); + let edges = vec![ + CompressedEdge { + vertices: (0, 1), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()), + (0.0, PI), + )), + }, + CompressedEdge { + vertices: (1, 0), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()), + (PI, 2.0 * PI), + )), + }, + CompressedEdge { + vertices: (0, 2), + curve: Curve::Line(Line(vertices[0], vertices[2])), + }, + CompressedEdge { + vertices: (2, 3), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed(translate), + (0.0, PI), + )), + }, + CompressedEdge { + vertices: (3, 2), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed(translate), + (PI, 2.0 * PI), + )), + }, + CompressedEdge { + vertices: (4, 5), + curve: Curve::PCurve(PCurve::new( + ParameterCurve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed(transform), + (0.0, PI), + )), + surface, + )), + }, + CompressedEdge { + vertices: (5, 4), + curve: Curve::PCurve(PCurve::new( + ParameterCurve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed(transform), + (PI, 2.0 * PI), + )), + surface, + )), + }, + ]; + let faces = vec![Face { + surface, + boundaries: vec![ + vec![ + CompressedEdgeIndex { + index: 1, + orientation: true, + }, + CompressedEdgeIndex { + index: 2, + orientation: true, + }, + CompressedEdgeIndex { + index: 4, + orientation: false, + }, + CompressedEdgeIndex { + index: 3, + orientation: false, + }, + CompressedEdgeIndex { + index: 2, + orientation: false, + }, + CompressedEdgeIndex { + index: 0, + orientation: true, + }, + ], + vec![ + CompressedEdgeIndex { + index: 6, + orientation: false, + }, + CompressedEdgeIndex { + index: 5, + orientation: false, + }, + ], + ], + orientation: true, + }]; + let mut shell = CompressedShell { + vertices, + edges, + faces, + }; + + assert!(Shell::extract(shell.clone()).is_err()); + split_closed_faces(&mut shell, 0.01, sp); + assert!(Shell::extract(shell.clone()).is_ok()); + + let CompressedShell { + ref vertices, + ref edges, + ref mut faces, + } = shell; + assert_eq!(vertices.len(), 6); + assert_eq!(edges.len(), 9); + assert_eq!(edges[7].vertices, (3, 5)); + let curve0 = &edges[7].curve; + let curve1 = Line(Point3::new(-1.0, 0.0, 1.0), Point3::new(-1.0, 0.0, 0.75)); + for i in 0..=10 { + let t = i as f64 / 10.0; + assert_near!(curve0.subs(t), curve1.subs(t)); + } + assert_eq!(edges[8].vertices, (4, 1)); + let curve0 = &edges[8].curve; + let curve1 = Line(Point3::new(-1.0, 0.0, 0.25), Point3::new(-1.0, 0.0, 0.0)); + for i in 0..=10 { + let t = i as f64 / 10.0; + assert_near!(curve0.subs(t), curve1.subs(t)); + } + assert_eq!(faces.len(), 2); + let i = faces + .iter_mut() + .position(|face| { + face.boundaries[0].contains(&CompressedEdgeIndex { + index: 2, + orientation: true, + }) + }) + .unwrap(); + if i == 1 { + faces.swap(0, 1); + } + let i = faces[0].boundaries[0] + .iter() + .position(|edge_index| { + *edge_index + == CompressedEdgeIndex { + index: 2, + orientation: true, + } + }) + .unwrap(); + faces[0].boundaries[0].rotate_left(i); + let i = faces[1].boundaries[0] + .iter() + .position(|edge_index| { + *edge_index + == CompressedEdgeIndex { + index: 3, + orientation: false, + } + }) + .unwrap(); + faces[1].boundaries[0].rotate_left(i); + + assert_eq!( + *faces, + vec![ + Face { + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 2, + orientation: true, + }, + CompressedEdgeIndex { + index: 4, + orientation: false, + }, + CompressedEdgeIndex { + index: 7, + orientation: true, + }, + CompressedEdgeIndex { + index: 5, + orientation: false, + }, + CompressedEdgeIndex { + index: 8, + orientation: true, + }, + CompressedEdgeIndex { + index: 1, + orientation: true, + } + ]], + surface, + orientation: true, + }, + Face { + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 3, + orientation: false, + }, + CompressedEdgeIndex { + index: 2, + orientation: false, + }, + CompressedEdgeIndex { + index: 0, + orientation: true, + }, + CompressedEdgeIndex { + index: 8, + orientation: false, + }, + CompressedEdgeIndex { + index: 6, + orientation: false, + }, + CompressedEdgeIndex { + index: 7, + orientation: false, + }, + ]], + surface, + orientation: true, + }, + ] + ); +} + +#[test] +fn test_split_closed_face_cylinder_with_rotated_hole() { + #[derive( + Clone, + Debug, + ParametricCurve, + BoundedCurve, + ParameterDivision1D, + Cut, + SearchNearestParameterD1, + )] + enum ParameterCurve { + Line(Line), + Arc(TrimmedCurve, Matrix3>>), + } + #[derive( + Clone, + Debug, + ParametricCurve, + BoundedCurve, + ParameterDivision1D, + Cut, + SearchNearestParameterD1, + )] + enum Curve { + Line(Line), + Arc(TrimmedCurve, Matrix4>>), + PCurve(PCurve), + } + impl From, Surface>> for Curve { + fn from(value: PCurve, Surface>) -> Self { + let (line, surface) = value.decompose(); + Self::PCurve(PCurve::new(ParameterCurve::Line(line), surface)) + } + } + type Surface = RevolutedCurve>; + + let surface = RevolutedCurve::by_revolution( + Line(Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 0.0, 0.0)), + Point3::origin(), + Vector3::unit_z(), + ); + + let vertices = vec![ + Point3::new(1.0, 0.0, 0.0), + Point3::new(-1.0, 0.0, 0.0), + Point3::new(1.0, 0.0, 1.0), + Point3::new(-1.0, 0.0, 1.0), + surface.subs(0.5, PI + 0.25), + surface.subs(0.5, PI - 0.25), + ]; + + let translate = Matrix4::from_translation(Vector3::unit_z()); + let transform = Matrix3::from_translation(Vector2::new(0.5, PI)) * Matrix3::from_scale(0.25); + let edges = vec![ + CompressedEdge { + vertices: (0, 1), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()), + (0.0, PI), + )), + }, + CompressedEdge { + vertices: (1, 0), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()), + (PI, 2.0 * PI), + )), + }, + CompressedEdge { + vertices: (0, 2), + curve: Curve::Line(Line(vertices[0], vertices[2])), + }, + CompressedEdge { + vertices: (2, 3), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed(translate), + (0.0, PI), + )), + }, + CompressedEdge { + vertices: (3, 2), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed(translate), + (PI, 2.0 * PI), + )), + }, + CompressedEdge { + vertices: (4, 5), + curve: Curve::PCurve(PCurve::new( + ParameterCurve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed(transform), + (0.5 * PI, 1.5 * PI), + )), + surface, + )), + }, + CompressedEdge { + vertices: (5, 4), + curve: Curve::PCurve(PCurve::new( + ParameterCurve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed(transform), + (1.5 * PI, 2.5 * PI), + )), + surface, + )), + }, + ]; + let faces = vec![Face { + surface, + boundaries: vec![ + vec![ + CompressedEdgeIndex { + index: 1, + orientation: true, + }, + CompressedEdgeIndex { + index: 2, + orientation: true, + }, + CompressedEdgeIndex { + index: 4, + orientation: false, + }, + CompressedEdgeIndex { + index: 3, + orientation: false, + }, + CompressedEdgeIndex { + index: 2, + orientation: false, + }, + CompressedEdgeIndex { + index: 0, + orientation: true, + }, + ], + vec![ + CompressedEdgeIndex { + index: 6, + orientation: false, + }, + CompressedEdgeIndex { + index: 5, + orientation: false, + }, + ], + ], + orientation: true, + }]; + let mut shell = CompressedShell { + vertices, + edges, + faces, + }; + + assert!(Shell::extract(shell.clone()).is_err()); + split_closed_faces(&mut shell, 0.01, sp); + assert!(Shell::extract(shell.clone()).is_ok()); + + let CompressedShell { + ref vertices, + ref edges, + ref mut faces, + } = shell; + assert_eq!(vertices.len(), 8); + assert_eq!(edges.len(), 11); + assert_eq!(edges[5].vertices, (4, 7)); + assert_eq!(edges[6].vertices, (5, 6)); + assert_eq!(edges[7].vertices, (6, 4)); + assert_eq!(edges[8].vertices, (7, 5)); + assert_eq!(edges[9].vertices, (3, 7)); + let curve0 = &edges[9].curve; + let curve1 = Line(Point3::new(-1.0, 0.0, 1.0), Point3::new(-1.0, 0.0, 0.75)); + for i in 0..=10 { + let t = i as f64 / 10.0; + assert_near!(curve0.subs(t), curve1.subs(t)); + } + assert_eq!(edges[10].vertices, (6, 1)); + let curve0 = &edges[10].curve; + let curve1 = Line(Point3::new(-1.0, 0.0, 0.25), Point3::new(-1.0, 0.0, 0.0)); + for i in 0..=10 { + let t = i as f64 / 10.0; + assert_near!(curve0.subs(t), curve1.subs(t)); + } + assert_eq!(faces.len(), 2); + let i = faces + .iter_mut() + .position(|face| { + face.boundaries[0].contains(&CompressedEdgeIndex { + index: 2, + orientation: true, + }) + }) + .unwrap(); + if i == 1 { + faces.swap(0, 1); + } + let i = faces[0].boundaries[0] + .iter() + .position(|edge_index| { + *edge_index + == CompressedEdgeIndex { + index: 2, + orientation: true, + } + }) + .unwrap(); + faces[0].boundaries[0].rotate_left(i); + let i = faces[1].boundaries[0] + .iter() + .position(|edge_index| { + *edge_index + == CompressedEdgeIndex { + index: 3, + orientation: false, + } + }) + .unwrap(); + faces[1].boundaries[0].rotate_left(i); + + assert_eq!( + *faces, + vec![ + Face { + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 2, + orientation: true, + }, + CompressedEdgeIndex { + index: 4, + orientation: false, + }, + CompressedEdgeIndex { + index: 9, + orientation: true, + }, + CompressedEdgeIndex { + index: 5, + orientation: false, + }, + CompressedEdgeIndex { + index: 7, + orientation: false, + }, + CompressedEdgeIndex { + index: 10, + orientation: true, + }, + CompressedEdgeIndex { + index: 1, + orientation: true, + } + ]], + surface, + orientation: true, + }, + Face { + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 3, + orientation: false, + }, + CompressedEdgeIndex { + index: 2, + orientation: false, + }, + CompressedEdgeIndex { + index: 0, + orientation: true, + }, + CompressedEdgeIndex { + index: 10, + orientation: false, + }, + CompressedEdgeIndex { + index: 6, + orientation: false, + }, + CompressedEdgeIndex { + index: 8, + orientation: false, + }, + CompressedEdgeIndex { + index: 9, + orientation: false, + }, + ]], + surface, + orientation: true, + }, + ] + ); +} + +#[test] +fn too_simple_cylinder() { + #[derive( + Clone, + Debug, + ParametricCurve, + BoundedCurve, + ParameterDivision1D, + Cut, + SearchNearestParameterD1, + )] + enum Curve { + Arc(TrimmedCurve, Matrix4>>), + PCurve(PCurve, Surface>), + } + impl From, Surface>> for Curve { + fn from(value: PCurve, Surface>) -> Self { Curve::PCurve(value) } + } + type Surface = RevolutedCurve>; + + let vertices = vec![Point3::new(1.0, 0.0, 0.0), Point3::new(1.0, 0.0, 1.0)]; + + let translation = Matrix4::from_translation(Vector3::unit_z()); + let circle0 = TrimmedCurve::new(Processor::new(UnitCircle::new()), (0.0, 2.0 * PI)); + let circle1 = TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed(translation), + (0.0, 2.0 * PI), + ); + let edges = vec![ + CompressedEdge { + vertices: (0, 0), + curve: Curve::Arc(circle0), + }, + CompressedEdge { + vertices: (1, 1), + curve: Curve::Arc(circle1), + }, + ]; + + let surface = RevolutedCurve::by_revolution( + Line(vertices[1], vertices[0]), + Point3::origin(), + Vector3::unit_z(), + ); + let faces = vec![CompressedFace { + boundaries: vec![ + vec![CompressedEdgeIndex { + index: 0, + orientation: true, + }], + vec![CompressedEdgeIndex { + index: 1, + orientation: false, + }], + ], + surface, + orientation: true, + }]; + + let mut shell = CompressedShell { + vertices, + edges, + faces, + }; + + assert!(Shell::extract(shell.clone()).is_err()); + split_closed_edges(&mut shell); + split_closed_faces(&mut shell, 0.01, sp); + assert!(Shell::extract(shell.clone()).is_ok()); + + let CompressedShell { + ref vertices, + ref edges, + ref mut faces, + } = shell; + + assert_eq!(vertices.len(), 4); + assert_near!(vertices[2], Point3::new(-1.0, 0.0, 0.0)); + assert_near!(vertices[3], Point3::new(-1.0, 0.0, 1.0)); + + assert_eq!(edges.len(), 6); + assert_eq!(edges[0].vertices, (0, 2)); + assert_eq!(edges[1].vertices, (1, 3)); + assert_eq!(edges[2].vertices, (2, 0)); + assert_eq!(edges[3].vertices, (3, 1)); + assert_eq!(edges[4].vertices, (0, 1)); + assert_eq!(edges[5].vertices, (2, 3)); + + assert_eq!(faces.len(), 2); + let i = faces + .iter_mut() + .position(|face| { + face.boundaries[0].contains(&CompressedEdgeIndex { + index: 0, + orientation: true, + }) + }) + .unwrap(); + if i == 1 { + faces.swap(0, 1); + } + let i = faces[0].boundaries[0] + .iter() + .position(|edge_index| { + *edge_index + == CompressedEdgeIndex { + index: 0, + orientation: true, + } + }) + .unwrap(); + faces[0].boundaries[0].rotate_left(i); + let i = faces[1].boundaries[0] + .iter() + .position(|edge_index| { + *edge_index + == CompressedEdgeIndex { + index: 2, + orientation: true, + } + }) + .unwrap(); + faces[1].boundaries[0].rotate_left(i); + + assert_eq!( + shell.faces, + vec![ + Face { + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 0, + orientation: true, + }, + CompressedEdgeIndex { + index: 5, + orientation: true, + }, + CompressedEdgeIndex { + index: 1, + orientation: false, + }, + CompressedEdgeIndex { + index: 4, + orientation: false, + }, + ]], + surface, + orientation: true, + }, + Face { + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 2, + orientation: true, + }, + CompressedEdgeIndex { + index: 4, + orientation: true, + }, + CompressedEdgeIndex { + index: 3, + orientation: false, + }, + CompressedEdgeIndex { + index: 5, + orientation: false, + }, + ]], + surface, + orientation: true, + }, + ] + ); +} + +#[test] +fn double_closed_boundary_cylinder() { + #[derive( + Clone, + Debug, + ParametricCurve, + BoundedCurve, + ParameterDivision1D, + Cut, + SearchNearestParameterD1, + )] + enum ParameterCurve { + Line(Line), + Arc(TrimmedCurve, Matrix3>>), + } + #[derive( + Clone, + Debug, + ParametricCurve, + BoundedCurve, + ParameterDivision1D, + Cut, + SearchNearestParameterD1, + )] + enum Curve { + Line(Line), + Arc(TrimmedCurve, Matrix4>>), + #[allow(clippy::enum_variant_names)] + PCurve(PCurve), + } + impl From, Surface>> for Curve { + fn from(value: PCurve, Surface>) -> Self { + let (line, surface) = value.decompose(); + Self::PCurve(PCurve::new(ParameterCurve::Line(line), surface)) + } + } + type Surface = RevolutedCurve>; + + let vertices = vec![ + Point3::new(1.0, 0.0, 0.0), + Point3::new(1.0, 0.0, 1.0), + Point3::new(1.0, 0.0, 0.25), + Point3::new(-1.0, 0.0, 0.25), + ]; + let surface = RevolutedCurve::by_revolution( + Line(vertices[1], vertices[0]), + Point3::origin(), + Vector3::unit_z(), + ); + let edges = vec![ + CompressedEdge { + vertices: (0, 0), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()), + (0.0, 2.0 * PI), + )), + }, + CompressedEdge { + vertices: (1, 1), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()) + .transformed(Matrix4::from_translation(Vector3::unit_z())), + (0.0, 2.0 * PI), + )), + }, + CompressedEdge { + vertices: (2, 2), + curve: Curve::PCurve(PCurve::new( + ParameterCurve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed( + Matrix3::from_translation(Vector2::new(0.5, 0.0)) + * Matrix3::from_scale(0.25), + ), + (0.0, 2.0 * PI), + )), + surface, + )), + }, + CompressedEdge { + vertices: (3, 3), + curve: Curve::PCurve(PCurve::new( + ParameterCurve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed( + Matrix3::from_translation(Vector2::new(0.5, PI)) + * Matrix3::from_scale(0.25), + ), + (0.0, 2.0 * PI), + )), + surface, + )), + }, + ]; + let faces = vec![CompressedFace { + boundaries: vec![ + vec![CompressedEdgeIndex { + index: 0, + orientation: true, + }], + vec![CompressedEdgeIndex { + index: 1, + orientation: false, + }], + vec![CompressedEdgeIndex { + index: 2, + orientation: false, + }], + vec![CompressedEdgeIndex { + index: 3, + orientation: false, + }], + ], + surface, + orientation: true, + }]; + let mut shell = CompressedShell { + vertices, + edges, + faces, + }; + + assert!(Shell::extract(shell.clone()).is_err()); + split_closed_edges(&mut shell); + split_closed_faces(&mut shell, 0.05, sp); + assert!(Shell::extract(shell.clone()).is_ok()); + + let CompressedShell { + ref vertices, + ref edges, + ref mut faces, + } = shell; + + assert_eq!(vertices.len(), 8); + assert_near!(vertices[0], Point3::new(1.0, 0.0, 0.0)); + assert_near!(vertices[1], Point3::new(1.0, 0.0, 1.0)); + assert_near!(vertices[2], Point3::new(1.0, 0.0, 0.25)); + assert_near!(vertices[3], Point3::new(-1.0, 0.0, 0.25)); + assert_near!(vertices[4], Point3::new(-1.0, 0.0, 0.0)); + assert_near!(vertices[5], Point3::new(-1.0, 0.0, 1.0)); + assert_near!(vertices[6], Point3::new(1.0, 0.0, 0.75)); + assert_near!(vertices[7], Point3::new(-1.0, 0.0, 0.75)); + + assert_eq!(edges.len(), 12); + assert_eq!(edges[0].vertices, (0, 4)); + assert_eq!(edges[1].vertices, (1, 5)); + assert_eq!(edges[2].vertices, (2, 6)); + assert_eq!(edges[3].vertices, (3, 7)); + assert_eq!(edges[4].vertices, (4, 0)); + assert_eq!(edges[5].vertices, (5, 1)); + assert_eq!(edges[6].vertices, (6, 2)); + assert_eq!(edges[7].vertices, (7, 3)); + assert_eq!(edges[8].vertices, (0, 2)); + assert_eq!(edges[9].vertices, (6, 1)); + assert_eq!(edges[10].vertices, (4, 3)); + assert_eq!(edges[11].vertices, (7, 5)); + + let i = faces + .iter_mut() + .position(|face| { + face.boundaries[0].contains(&CompressedEdgeIndex { + index: 0, + orientation: true, + }) + }) + .unwrap(); + if i == 1 { + faces.swap(0, 1); + } + let i = faces[0].boundaries[0] + .iter() + .position(|edge_index| { + *edge_index + == CompressedEdgeIndex { + index: 0, + orientation: true, + } + }) + .unwrap(); + faces[0].boundaries[0].rotate_left(i); + let i = faces[1].boundaries[0] + .iter() + .position(|edge_index| { + *edge_index + == CompressedEdgeIndex { + index: 4, + orientation: true, + } + }) + .unwrap(); + faces[1].boundaries[0].rotate_left(i); + + assert_eq!( + *faces, + vec![ + CompressedFace { + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 0, + orientation: true, + }, + CompressedEdgeIndex { + index: 10, + orientation: true, + }, + CompressedEdgeIndex { + index: 7, + orientation: false, + }, + CompressedEdgeIndex { + index: 11, + orientation: true, + }, + CompressedEdgeIndex { + index: 1, + orientation: false, + }, + CompressedEdgeIndex { + index: 9, + orientation: false, + }, + CompressedEdgeIndex { + index: 2, + orientation: false, + }, + CompressedEdgeIndex { + index: 8, + orientation: false, + }, + ]], + surface, + orientation: true, + }, + CompressedFace { + boundaries: vec![vec![ + CompressedEdgeIndex { + index: 4, + orientation: true, + }, + CompressedEdgeIndex { + index: 8, + orientation: true, + }, + CompressedEdgeIndex { + index: 6, + orientation: false, + }, + CompressedEdgeIndex { + index: 9, + orientation: true, + }, + CompressedEdgeIndex { + index: 5, + orientation: false, + }, + CompressedEdgeIndex { + index: 11, + orientation: false, + }, + CompressedEdgeIndex { + index: 3, + orientation: false, + }, + CompressedEdgeIndex { + index: 10, + orientation: false, + }, + ]], + surface, + orientation: true, + }, + ] + ); +} + +#[test] +fn many_closed_boundary_cylinder() { + #[derive( + Clone, + Debug, + ParametricCurve, + BoundedCurve, + ParameterDivision1D, + Cut, + SearchNearestParameterD1, + )] + enum ParameterCurve { + Line(Line), + Arc(TrimmedCurve, Matrix3>>), + } + #[derive( + Clone, + Debug, + ParametricCurve, + BoundedCurve, + ParameterDivision1D, + Cut, + SearchNearestParameterD1, + )] + enum Curve { + Line(Line), + Arc(TrimmedCurve, Matrix4>>), + PCurve(PCurve), + } + impl From, Surface>> for Curve { + fn from(value: PCurve, Surface>) -> Self { + let (line, surface) = value.decompose(); + Self::PCurve(PCurve::new(ParameterCurve::Line(line), surface)) + } + } + type Surface = RevolutedCurve>; + + const NUM_OF_CIRCLES: usize = 10; + + let vertices = { + let mut vertices = vec![Point3::new(1.0, 0.0, 0.0), Point3::new(1.0, 0.0, 1.0)]; + vertices.extend((0..NUM_OF_CIRCLES).map(|i| { + let t = 2.0 * PI * i as f64 / NUM_OF_CIRCLES as f64; + Point3::new(f64::cos(t), f64::sin(t), 0.4) + })); + vertices + }; + let surface = RevolutedCurve::by_revolution( + Line(vertices[1], vertices[0]), + Point3::origin(), + Vector3::unit_z(), + ); + let edges = { + let mut edges = vec![ + CompressedEdge { + vertices: (0, 0), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()), + (0.0, 2.0 * PI), + )), + }, + CompressedEdge { + vertices: (1, 1), + curve: Curve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()) + .transformed(Matrix4::from_translation(Vector3::unit_z())), + (0.0, 2.0 * PI), + )), + }, + ]; + edges.extend((0..NUM_OF_CIRCLES).map(|i| { + let t = 2.0 * PI * i as f64 / NUM_OF_CIRCLES as f64; + CompressedEdge { + vertices: (2 + i, 2 + i), + curve: Curve::PCurve(PCurve::new( + ParameterCurve::Arc(TrimmedCurve::new( + Processor::new(UnitCircle::new()).transformed( + Matrix3::from_translation(Vector2::new(0.5, t)) + * Matrix3::from_scale(0.1), + ), + (0.0, 2.0 * PI), + )), + surface, + )), + } + })); + edges + }; + let mut boundaries = vec![vec![CompressedEdgeIndex { + index: 0, + orientation: true, + }]]; + boundaries.extend((0..=NUM_OF_CIRCLES).map(|i| { + vec![CompressedEdgeIndex { + index: 1 + i, + orientation: false, + }] + })); + let faces = vec![CompressedFace { + boundaries, + surface, + orientation: true, + }]; + let mut shell = CompressedShell { + vertices, + edges, + faces, + }; + + assert!(Shell::extract(shell.clone()).is_err()); + split_closed_edges(&mut shell); + split_closed_faces(&mut shell, 0.05, sp); + assert!(Shell::extract(shell.clone()).is_ok()); + + let CompressedShell { + ref vertices, + ref edges, + ref mut faces, + } = shell; + + assert_eq!(vertices.len(), (2 + NUM_OF_CIRCLES) * 2); + assert_near!(vertices[0], Point3::new(1.0, 0.0, 0.0)); + assert_near!(vertices[1], Point3::new(1.0, 0.0, 1.0)); + assert_near!(vertices[NUM_OF_CIRCLES + 2], Point3::new(-1.0, 0.0, 0.0)); + assert_near!(vertices[NUM_OF_CIRCLES + 3], Point3::new(-1.0, 0.0, 1.0)); + (0..NUM_OF_CIRCLES).for_each(|i| { + let t = 2.0 * PI * i as f64 / NUM_OF_CIRCLES as f64; + assert_near!(vertices[2 + i], Point3::new(f64::cos(t), f64::sin(t), 0.4)); + assert_near!( + vertices[4 + i + NUM_OF_CIRCLES], + Point3::new(f64::cos(t), f64::sin(t), 0.6) + ); + }); + + assert_eq!(edges.len(), 8 + NUM_OF_CIRCLES * 2); + (0..NUM_OF_CIRCLES + 2).for_each(|i| { + let j = i + 2 + NUM_OF_CIRCLES; + assert_eq!(edges[i].vertices, (i, j)); + assert_eq!(edges[j].vertices, (j, i)); + }); + let i = 4 + NUM_OF_CIRCLES * 2; + assert_eq!(edges[i].vertices, (0, 2)); + assert_eq!(edges[i + 1].vertices, (4 + NUM_OF_CIRCLES, 1)); + assert_eq!( + edges[i + 2].vertices, + (2 + NUM_OF_CIRCLES, 2 + NUM_OF_CIRCLES / 2) + ); + assert_eq!( + edges[i + 3].vertices, + (4 + NUM_OF_CIRCLES * 3 / 2, 3 + NUM_OF_CIRCLES) + ); + + let i = faces + .iter_mut() + .position(|face| { + face.boundaries[0].contains(&CompressedEdgeIndex { + index: 0, + orientation: true, + }) + }) + .unwrap(); + if i == 1 { + faces.swap(0, 1); + } + + assert_eq!(faces[0].boundaries.len(), NUM_OF_CIRCLES / 2); + assert_eq!(faces[1].boundaries.len(), NUM_OF_CIRCLES / 2); + + (1..NUM_OF_CIRCLES / 2).for_each(|i| { + (0..=1).for_each(|fid| { + let eid = faces[fid].boundaries[i][0].index; + let p = vertices[edges[eid].vertices.0]; + assert!(f64::signum(p.y) == f64::signum(f64::powi(-1.0, fid as i32))); + }); + }); +} + +fn sp(surface: &S, p: Point3, hint: Option<(f64, f64)>) -> Option<(f64, f64)> +where S: SearchParameter { + surface.search_parameter(p, hint, 10) +} + +#[test] +#[cfg(feature = "step-test")] +fn step_import() { + use truck_stepio::r#in::*; + const STEP_DIRECTORY: &str = concat!(env!("CARGO_MANIFEST_DIR"), "/../resources/step/"); + const STEP_FILES: &[&str] = &[ + "occt-cylinder.step", + "occt-cone.step", + "abc-0006.step", + "abc-0008.step", + "abc-0035.step", + ]; + + STEP_FILES.iter().for_each(|file_name| { + println!("{file_name}"); + let path = [STEP_DIRECTORY, file_name].concat(); + let step_string = std::fs::read_to_string(path).unwrap(); + let table = Table::from_step(&step_string).unwrap(); + table.shell.values().cloned().for_each(|step_shell| { + let mut cshell = table.to_compressed_shell(&step_shell).unwrap(); + cshell.robust_split_closed_edges_and_faces(0.05); + truck_topology::Shell::extract(cshell).unwrap(); + }); + }); +} diff --git a/vendor/truck-shapeops/src/lib.rs b/vendor/truck-shapeops/src/lib.rs new file mode 100644 index 00000000..dabc8657 --- /dev/null +++ b/vendor/truck-shapeops/src/lib.rs @@ -0,0 +1,37 @@ +//! Crate for operationg shapes. Provides boolean operations to Solid, and shape healing for importing shapes from other CAD systems. + +// SolveCraft: warnings are not errors in the vendored copy (newer compilers flag more). +#![allow(dead_code)] +#![warn(clippy::all, rust_2018_idioms)] +#![warn( + missing_docs, + missing_debug_implementations, + trivial_casts, + trivial_numeric_casts, + unsafe_code, + unstable_features, + unused_import_braces, + unused_qualifications +)] + +/// SolveCraft: report where an operation gives up (set `SHAPEOPS_TRACE`). +macro_rules! tr { + ($e:expr, $tag:expr) => { + match $e { + Some(v) => v, + None => { + if std::env::var_os("SHAPEOPS_TRACE").is_some() { + eprintln!("shapeops: gave up at {} ({}:{})", $tag, file!(), line!()); + } + return None; + } + } + }; +} +pub(crate) use tr; + +mod healing; +pub use healing::{RobustSplitClosedEdgesAndFaces, SplitClosedEdgesAndFaces}; +mod transversal; +pub use transversal::{and, or, ShapeOpsCurve, ShapeOpsSurface}; +mod alternative; diff --git a/vendor/truck-shapeops/src/transversal/divide_face/mod.rs b/vendor/truck-shapeops/src/transversal/divide_face/mod.rs new file mode 100644 index 00000000..e72cef0e --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/divide_face/mod.rs @@ -0,0 +1,140 @@ +#![allow(clippy::many_single_char_names)] + +use super::faces_classification::FacesClassification; +use super::loops_store::*; +use rustc_hash::FxHashMap as HashMap; +use std::ops::Deref; +use truck_meshalgo::prelude::*; +use truck_topology::*; + +fn create_parameter_boundary( + face: &Face, + wire: &Wire, + polys: &mut HashMap, PolylineCurve

>, + tol: f64, +) -> Option> +where + P: Copy, + C: BoundedCurve + ParameterDivision1D, + S: Clone + SearchParameter + SearchNearestParameter, +{ + let surface = face.surface(); + let pt = wire.front_vertex().unwrap().point(); + // SolveCraft: intersection-curve points lie on the surface only within the tolerance; fall + // back to the nearest parameter when the exact search fails. + let p: Point2 = surface + .search_parameter(pt, None, 100) + .or_else(|| surface.search_nearest_parameter(pt, None, 100))? + .into(); + let vec = wire.edge_iter().try_fold(vec![p], |mut vec, edge| { + let poly = polys.entry(edge.id()).or_insert_with(|| { + let curve = edge.curve(); + let div = curve.parameter_division(curve.range_tuple(), tol).1; + PolylineCurve(div) + }); + let mut p = *vec.last().unwrap(); + let closure = |q: &P| -> Option { + p = surface + .search_parameter(*q, Some(p.into()), 100) + .or_else(|| surface.search_nearest_parameter(*q, Some(p.into()), 100))? + .into(); + Some(p) + }; + let add: Option> = match edge.orientation() { + true => poly.iter().skip(1).map(closure).collect(), + false => poly.iter().rev().skip(1).map(closure).collect(), + }; + vec.append(&mut add?); + Some(vec) + })?; + Some(PolylineCurve(vec)) +} + +#[derive(Clone, Debug)] +struct WireChunk<'a, C> { + poly: PolylineCurve, + wire: &'a BoundaryWire, +} + +type FaceWithShapesOpStatus = (Face, ShapesOpStatus); +fn divide_one_face( + face: &Face, + loops: &Loops, + tol: f64, +) -> Option>> +where + C: BoundedCurve + ParameterDivision1D, + S: Clone + SearchParameter + SearchNearestParameter, +{ + let (mut pre_faces, mut negative_wires) = (Vec::new(), Vec::new()); + let mut map = HashMap::default(); + loops.iter().try_for_each(|wire| { + let poly = crate::tr!(create_parameter_boundary(face, wire, &mut map, tol), "parameter boundary"); + match poly.area() > 0.0 { + true => pre_faces.push(vec![WireChunk { poly, wire }]), + false => negative_wires.push(WireChunk { poly, wire }), + } + Some(()) + })?; + negative_wires.into_iter().try_for_each(|chunk| { + let pt = chunk.poly.front(); + let op = crate::tr!(pre_faces.iter_mut().find(|face| face[0].poly.include(pt)), "loop not inside any face piece"); + op.push(chunk); + Some(()) + })?; + let vec: Vec<_> = pre_faces + .into_iter() + .map(|pre_face| { + let surface = face.surface(); + let op = pre_face + .iter() + .find(|chunk| chunk.wire.status() != ShapesOpStatus::Unknown); + let status = match op { + Some(chunk) => chunk.wire.status(), + None => ShapesOpStatus::Unknown, + }; + let wires: Vec> = pre_face + .into_iter() + .map(|chunk| chunk.wire.deref().clone()) + .collect(); + let mut new_face = Face::debug_new(wires, surface); + if !face.orientation() { + new_face.invert(); + } + (new_face, status) + }) + .collect(); + Some(vec) +} + +pub fn divide_faces( + shell: &Shell, + loops_store: &LoopsStore, + tol: f64, +) -> Option> +where + C: BoundedCurve + ParameterDivision1D, + S: Clone + SearchParameter + SearchNearestParameter, +{ + let mut res = FacesClassification::::default(); + shell + .iter() + .zip(loops_store) + .try_for_each(|(face, loops)| { + if loops + .iter() + .all(|wire| wire.status() == ShapesOpStatus::Unknown) + { + res.push(face.clone(), ShapesOpStatus::Unknown); + } else { + let vec = crate::tr!(divide_one_face(face, loops, tol), "divide one face"); + vec.into_iter() + .for_each(|(face, status)| res.push(face, status)); + } + Some(()) + })?; + Some(res) +} + +#[cfg(test)] +mod tests; diff --git a/vendor/truck-shapeops/src/transversal/divide_face/tests.rs b/vendor/truck-shapeops/src/transversal/divide_face/tests.rs new file mode 100644 index 00000000..7856330d --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/divide_face/tests.rs @@ -0,0 +1,226 @@ +use super::*; +use shell::ShellCondition; +use truck_geometry::prelude::*; +use truck_topology::Vertex; +const TOL: f64 = 0.05; + +fn line(v0: &Vertex, v1: &Vertex) -> Edge> { + let curve = BSplineCurve::new(KnotVec::bezier_knot(1), vec![v0.point(), v1.point()]); + Edge::new(v0, v1, curve) +} + +fn parabola( + v0: &Vertex, + v1: &Vertex, + pt: Point3, +) -> Edge> { + let curve = BSplineCurve::new(KnotVec::bezier_knot(2), vec![v0.point(), pt, v1.point()]); + Edge::new(v0, v1, curve) +} + +#[test] +fn divide_plane_test() { + let v = Vertex::news([ + Point3::new(0.0, 0.0, 0.0), + Point3::new(0.0, 4.0, 0.0), + Point3::new(-1.0, 1.0, 0.0), + Point3::new(-1.0, 3.0, 0.0), + Point3::new(1.0, 1.0, 0.0), + Point3::new(1.0, 3.0, 0.0), + ]); + let edge = vec![ + parabola(&v[0], &v[1], Point3::new(-4.0, 2.0, 0.0)), + parabola(&v[0], &v[1], Point3::new(4.0, 2.0, 0.0)), + line(&v[0], &v[1]), + parabola(&v[2], &v[3], Point3::new(-3.0, 2.0, 0.0)), + parabola(&v[2], &v[3], Point3::new(-1.0, 2.0, 0.0)), + parabola(&v[4], &v[5], Point3::new(1.0, 2.0, 0.0)), + parabola(&v[4], &v[5], Point3::new(3.0, 2.0, 0.0)), + ]; + let wire: Vec> = vec![ + vec![edge[1].clone(), edge[0].inverse()].into(), + vec![edge[2].clone(), edge[0].inverse()].into(), + vec![edge[1].clone(), edge[2].inverse()].into(), + vec![edge[3].clone(), edge[4].inverse()].into(), + vec![edge[5].clone(), edge[6].inverse()].into(), + ]; + let face = Face::new( + vec![wire[0].clone(), wire[3].clone(), wire[4].clone()], + Plane::new( + Point3::origin(), + Point3::new(1.0, 0.0, 0.0), + Point3::new(0.0, 1.0, 0.0), + ), + ); + let loops: Loops<_, _> = vec![ + BoundaryWire::new(wire[1].clone(), ShapesOpStatus::Or), + BoundaryWire::new(wire[2].clone(), ShapesOpStatus::And), + BoundaryWire::new(wire[3].clone(), ShapesOpStatus::Unknown), + BoundaryWire::new(wire[4].clone(), ShapesOpStatus::Unknown), + ] + .into_iter() + .collect(); + let res = divide_one_face(&face, &loops, 0.01).unwrap(); + assert_eq!(res.len(), 2); + let (mut or, mut and) = (true, true); + for (face, status) in res { + let bdd = face.absolute_boundaries(); + match status { + ShapesOpStatus::Or => { + assert_eq!(bdd.len(), 2); + assert!(bdd[0] == wire[1] || bdd[0] == wire[3]); + assert!(bdd[1] == wire[1] || bdd[1] == wire[3]); + assert_ne!(bdd[0], bdd[1]); + assert!(or); + or = false; + } + ShapesOpStatus::And => { + assert_eq!(bdd.len(), 2); + assert!(bdd[0] == wire[2] || bdd[0] == wire[4]); + assert!(bdd[1] == wire[2] || bdd[1] == wire[4]); + assert_ne!(bdd[0], bdd[1]); + assert!(and); + and = false; + } + _ => panic!("There must be no unknown!"), + } + } +} + +type AlternativeIntersection = crate::alternative::Alternative< + NurbsCurve, + IntersectionCurve, AlternativeSurface>, +>; +type AlternativeSurface = crate::alternative::Alternative, Plane>; + +crate::impl_from!( + NurbsCurve, + IntersectionCurve, AlternativeSurface> +); +crate::impl_from!(BSplineSurface, Plane); + +fn parabola_surfaces() -> (AlternativeSurface, AlternativeSurface) { + // define surfaces + #[rustfmt::skip] + let ctrl0 = vec![ + vec![Point3::new(-1.0, -1.0, 3.0), Point3::new(-1.0, 0.0, -1.0), Point3::new(-1.0, 1.0, 3.0)], + vec![Point3::new(0.0, -1.0, -1.0), Point3::new(0.0, 0.0, -5.0), Point3::new(0.0, 1.0, -1.0)], + vec![Point3::new(1.0, -1.0, 3.0), Point3::new(1.0, 0.0, -1.0), Point3::new(1.0, 1.0, 3.0)], + ]; + #[rustfmt::skip] + let ctrl1 = vec![ + vec![Point3::new(-1.0, -1.0, -3.0), Point3::new(-1.0, 0.0, 1.0), Point3::new(-1.0, 1.0, -3.0)], + vec![Point3::new(0.0, -1.0, 1.0), Point3::new(0.0, 0.0, 5.0), Point3::new(0.0, 1.0, 1.0)], + vec![Point3::new(1.0, -1.0, -3.0), Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 1.0, -3.0)], + ]; + ( + BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl0).into(), + BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl1).into(), + ) +} + +#[test] +fn independent_intersection() { + // prepare geoetries + let arc00: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, 1.0, 0.0, 0.0), + Vector4::new(-1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc01: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(-1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, -1.0, 0.0, 0.0), + Vector4::new(1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc10: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, -1.0, 1.0), + Vector4::new(0.0, 1.0, 0.0, 0.0), + Vector4::new(-1.0, 0.0, -1.0, 1.0), + ], + )) + .into(); + let arc11: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(-1.0, 0.0, -1.0, 1.0), + Vector4::new(0.0, -1.0, 0.0, 0.0), + Vector4::new(1.0, 0.0, -1.0, 1.0), + ], + )) + .into(); + let (surface0, surface1) = parabola_surfaces(); + let plane0: AlternativeSurface = Plane::new( + Point3::new(0.0, 0.0, 1.0), + Point3::new(1.0, 0.0, 1.0), + Point3::new(0.0, 1.0, 1.0), + ) + .into(); + let plane1: AlternativeSurface = Plane::new( + Point3::new(0.0, 0.0, -1.0), + Point3::new(1.0, 0.0, -1.0), + Point3::new(0.0, 1.0, -1.0), + ) + .into(); + + // prepare topologies + let v00 = Vertex::new(Point3::new(1.0, 0.0, 1.0)); + let v01 = Vertex::new(Point3::new(-1.0, 0.0, 1.0)); + let v10 = Vertex::new(Point3::new(1.0, 0.0, -1.0)); + let v11 = Vertex::new(Point3::new(-1.0, 0.0, -1.0)); + let wire0: Wire<_, _> = vec![Edge::new(&v00, &v01, arc00), Edge::new(&v01, &v00, arc01)].into(); + let wire1: Wire<_, _> = vec![Edge::new(&v10, &v11, arc10), Edge::new(&v11, &v10, arc11)].into(); + let shell0: Shell<_, _, _> = vec![ + Face::new(vec![wire0.clone()], plane0), + Face::new(vec![wire0], surface0).inverse(), + ] + .into(); + assert_eq!(shell0.shell_condition(), ShellCondition::Closed); + let shell1: Shell<_, _, _> = vec![ + Face::new(vec![wire1.clone()], plane1).inverse(), + Face::new(vec![wire1], surface1), + ] + .into(); + assert_eq!(shell1.shell_condition(), ShellCondition::Closed); + let poly_shell0 = shell0.triangulation(TOL); + let poly_shell1 = shell1.triangulation(TOL); + + let LoopsStoreQuadruple { + geom_loops_store0: loops_store0, + geom_loops_store1: loops_store1, + .. + } = create_loops_stores(&shell0, &poly_shell0, &shell1, &poly_shell1, TOL).unwrap(); + let [and0, or0, unknown0] = divide_faces(&shell0, &loops_store0, TOL) + .unwrap() + .and_or_unknown(); + let [and1, or1, unknown1] = divide_faces(&shell1, &loops_store1, TOL) + .unwrap() + .and_or_unknown(); + assert_eq!(and0.len(), 1); + assert_eq!(or0.len(), 1); + assert_eq!(unknown0.len(), 1); + assert_eq!(and1.len(), 1); + assert_eq!(or1.len(), 1); + assert_eq!(unknown1.len(), 1); + + match unknown0[0].surface() { + AlternativeSurface::FirstType(_) => panic!("This face must plane!"), + AlternativeSurface::SecondType(_) => {} + } + match unknown1[0].surface() { + AlternativeSurface::FirstType(_) => panic!("This face must plane!"), + AlternativeSurface::SecondType(_) => {} + } + + let and_shell: Shell<_, _, _> = vec![and0[0].clone(), and1[0].clone()].into(); + assert_eq!(and_shell.shell_condition(), ShellCondition::Closed); +} diff --git a/vendor/truck-shapeops/src/transversal/faces_classification/mod.rs b/vendor/truck-shapeops/src/transversal/faces_classification/mod.rs new file mode 100644 index 00000000..aeb674aa --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/faces_classification/mod.rs @@ -0,0 +1,67 @@ +use super::loops_store::ShapesOpStatus; +use rustc_hash::FxHashMap as HashMap; +use truck_topology::*; + +#[derive(Clone, Debug)] +pub struct FacesClassification { + shell: Shell, + status: HashMap, ShapesOpStatus>, +} + +impl Default for FacesClassification { + fn default() -> Self { + Self { + shell: Default::default(), + status: HashMap::default(), + } + } +} + +impl FacesClassification { + pub fn push(&mut self, face: Face, status: ShapesOpStatus) { + self.status.insert(face.id(), status); + self.shell.push(face); + } + + pub fn and_or_unknown(&self) -> [Shell; 3] { + let [mut and, mut or, mut unknown] = <[Shell; 3]>::default(); + for face in &self.shell { + match self.status.get(&face.id()).unwrap() { + ShapesOpStatus::And => and.push(face.clone()), + ShapesOpStatus::Or => or.push(face.clone()), + ShapesOpStatus::Unknown => unknown.push(face.clone()), + } + } + [and, or, unknown] + } + + pub fn integrate_by_component(&mut self) { + let [and, or, unknown] = self.and_or_unknown(); + let and_boundary = and.extract_boundaries(); + let or_boundary = or.extract_boundaries(); + let components = unknown.connected_components(); + for comp in components { + let boundary = comp.extract_boundaries(); + if and_boundary + .iter() + .flatten() + .any(|edge| edge.id() == boundary[0][0].id()) + { + comp.iter().for_each(|face| { + *self.status.get_mut(&face.id()).unwrap() = ShapesOpStatus::And; + }) + } else if or_boundary + .iter() + .flatten() + .any(|edge| edge.id() == boundary[0][0].id()) + { + comp.iter().for_each(|face| { + *self.status.get_mut(&face.id()).unwrap() = ShapesOpStatus::Or; + }) + } + } + } +} + +#[cfg(test)] +mod tests; diff --git a/vendor/truck-shapeops/src/transversal/faces_classification/tests.rs b/vendor/truck-shapeops/src/transversal/faces_classification/tests.rs new file mode 100644 index 00000000..a839602d --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/faces_classification/tests.rs @@ -0,0 +1,131 @@ +use super::{ + super::{divide_face, loops_store}, + *, +}; +use truck_geometry::prelude::*; +use truck_meshalgo::prelude::*; +use truck_topology::{shell::ShellCondition, Vertex}; +const TOL: f64 = 0.05; + +type AlternativeIntersection = crate::alternative::Alternative< + NurbsCurve, + IntersectionCurve, AlternativeSurface>, +>; +type AlternativeSurface = crate::alternative::Alternative, Plane>; + +fn parabola_surfaces() -> (AlternativeSurface, AlternativeSurface) { + // define surfaces + #[rustfmt::skip] + let ctrl0 = vec![ + vec![Point3::new(-1.0, -1.0, 3.0), Point3::new(-1.0, 0.0, -1.0), Point3::new(-1.0, 1.0, 3.0)], + vec![Point3::new(0.0, -1.0, -1.0), Point3::new(0.0, 0.0, -5.0), Point3::new(0.0, 1.0, -1.0)], + vec![Point3::new(1.0, -1.0, 3.0), Point3::new(1.0, 0.0, -1.0), Point3::new(1.0, 1.0, 3.0)], + ]; + #[rustfmt::skip] + let ctrl1 = vec![ + vec![Point3::new(-1.0, -1.0, -3.0), Point3::new(-1.0, 0.0, 1.0), Point3::new(-1.0, 1.0, -3.0)], + vec![Point3::new(0.0, -1.0, 1.0), Point3::new(0.0, 0.0, 5.0), Point3::new(0.0, 1.0, 1.0)], + vec![Point3::new(1.0, -1.0, -3.0), Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 1.0, -3.0)], + ]; + ( + BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl0).into(), + BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl1).into(), + ) +} + +#[test] +fn independent_intersection() { + // prepare geoetries + let arc00: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, 1.0, 0.0, 0.0), + Vector4::new(-1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc01: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(-1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, -1.0, 0.0, 0.0), + Vector4::new(1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc10: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, -1.0, 1.0), + Vector4::new(0.0, 1.0, 0.0, 0.0), + Vector4::new(-1.0, 0.0, -1.0, 1.0), + ], + )) + .into(); + let arc11: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(-1.0, 0.0, -1.0, 1.0), + Vector4::new(0.0, -1.0, 0.0, 0.0), + Vector4::new(1.0, 0.0, -1.0, 1.0), + ], + )) + .into(); + let (surface0, surface1) = parabola_surfaces(); + let plane0: AlternativeSurface = Plane::new( + Point3::new(0.0, 0.0, 1.0), + Point3::new(1.0, 0.0, 1.0), + Point3::new(0.0, 1.0, 1.0), + ) + .into(); + let plane1: AlternativeSurface = Plane::new( + Point3::new(0.0, 0.0, -1.0), + Point3::new(1.0, 0.0, -1.0), + Point3::new(0.0, 1.0, -1.0), + ) + .into(); + + // prepare topologies + let v00 = Vertex::new(Point3::new(1.0, 0.0, 1.0)); + let v01 = Vertex::new(Point3::new(-1.0, 0.0, 1.0)); + let v10 = Vertex::new(Point3::new(1.0, 0.0, -1.0)); + let v11 = Vertex::new(Point3::new(-1.0, 0.0, -1.0)); + let wire0: Wire<_, _> = vec![Edge::new(&v00, &v01, arc00), Edge::new(&v01, &v00, arc01)].into(); + let wire1: Wire<_, _> = vec![Edge::new(&v10, &v11, arc10), Edge::new(&v11, &v10, arc11)].into(); + let shell0: Shell<_, _, _> = vec![ + Face::new(vec![wire0.clone()], plane0), + Face::new(vec![wire0], surface0).inverse(), + ] + .into(); + assert_eq!(shell0.shell_condition(), ShellCondition::Closed); + let shell1: Shell<_, _, _> = vec![ + Face::new(vec![wire1.clone()], plane1).inverse(), + Face::new(vec![wire1], surface1), + ] + .into(); + assert_eq!(shell1.shell_condition(), ShellCondition::Closed); + let poly_shell0 = shell0.triangulation(TOL); + let poly_shell1 = shell1.triangulation(TOL); + + let loops_store::LoopsStoreQuadruple { + geom_loops_store0: loops_store0, + geom_loops_store1: loops_store1, + .. + } = loops_store::create_loops_stores(&shell0, &poly_shell0, &shell1, &poly_shell1, TOL) + .unwrap(); + let mut cls0 = divide_face::divide_faces(&shell0, &loops_store0, TOL).unwrap(); + cls0.integrate_by_component(); + let mut cls1 = divide_face::divide_faces(&shell1, &loops_store1, TOL).unwrap(); + cls1.integrate_by_component(); + + let [mut and, mut or, _] = cls0.and_or_unknown(); + let [and1, or1, _] = cls1.and_or_unknown(); + and.extend(and1); + or.extend(or1); + + assert_eq!(and.len(), 2); + assert_eq!(or.len(), 4); + assert_eq!(and.shell_condition(), ShellCondition::Closed); + assert_eq!(or.shell_condition(), ShellCondition::Closed); +} diff --git a/vendor/truck-shapeops/src/transversal/integrate/mod.rs b/vendor/truck-shapeops/src/transversal/integrate/mod.rs new file mode 100644 index 00000000..8e5d1248 --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/integrate/mod.rs @@ -0,0 +1,151 @@ +use super::*; +use truck_geometry::prelude::IntersectionCurve; +use truck_meshalgo::prelude::*; +use truck_topology::*; + +/// Only solids consisting of faces whose surface is implemented this trait can be used for set operations. +pub trait ShapeOpsSurface: + ParametricSurface3D + + ParameterDivision2D + + SearchParameter + + SearchNearestParameter + + Invertible + + Send + + Sync { +} +impl ShapeOpsSurface for S where S: ParametricSurface3D + + ParameterDivision2D + + SearchParameter + + SearchNearestParameter + + Invertible + + Send + + Sync +{ +} + +/// Only solids consisting of edges whose curve is implemented this trait can be used for set operations. +pub trait ShapeOpsCurve: + ParametricCurve3D + + ParameterDivision1D + + Cut + + Invertible + + From, S>> + + SearchParameter + + SearchNearestParameter + + Send + + Sync { +} +impl ShapeOpsCurve for C where C: ParametricCurve3D + + ParameterDivision1D + + Cut + + Invertible + + From, S>> + + SearchParameter + + SearchNearestParameter + + Send + + Sync +{ +} + +fn process_one_pair_of_shells, S: ShapeOpsSurface>( + shell0: &Shell, + shell1: &Shell, + tol: f64, +) -> Option<[Shell; 2]> { + nonpositive_tolerance!(tol); + // SolveCraft: robust meshing (boundaries from earlier operations lie on surfaces only + // within tolerance). + let poly_shell0 = shell0.robust_triangulation(tol); + let poly_shell1 = shell1.robust_triangulation(tol); + let loops_store::LoopsStoreQuadruple { + geom_loops_store0: loops_store0, + geom_loops_store1: loops_store1, + .. + } = crate::tr!(loops_store::create_loops_stores(shell0, &poly_shell0, shell1, &poly_shell1, tol), "loops stores"); + let mut cls0 = crate::tr!(divide_face::divide_faces(shell0, &loops_store0, tol), "divide faces 0"); + cls0.integrate_by_component(); + let mut cls1 = crate::tr!(divide_face::divide_faces(shell1, &loops_store1, tol), "divide faces 1"); + cls1.integrate_by_component(); + let [mut and0, mut or0, unknown0] = cls0.and_or_unknown(); + unknown0.into_iter().try_for_each(|face| { + let pt = face.boundaries()[0].vertex_iter().next().unwrap().point(); + let dir = hash::take_one_unit(pt); + let count = poly_shell1.iter().try_fold(0, |count, face| { + let poly = face.surface()?; + Some(count + poly.signed_crossing_faces(pt, dir)) + })?; + if count >= 1 { + and0.push(face); + } else { + or0.push(face); + } + Some(()) + })?; + let [mut and1, mut or1, unknown1] = cls1.and_or_unknown(); + unknown1.into_iter().try_for_each(|face| { + let pt = face.boundaries()[0].vertex_iter().next().unwrap().point(); + let dir = hash::take_one_unit(pt); + let count = poly_shell0.iter().try_fold(0, |count, face| { + let poly = face.surface()?; + Some(count + poly.signed_crossing_faces(pt, dir)) + })?; + if count >= 1 { + and1.push(face); + } else { + or1.push(face); + } + Some(()) + })?; + and0.append(&mut and1); + or0.append(&mut or1); + Some([and0, or0]) +} + +/// AND operation between two solids. +pub fn and, S: ShapeOpsSurface>( + solid0: &Solid, + solid1: &Solid, + tol: f64, +) -> Option> { + let mut iter0 = solid0.boundaries().iter(); + let mut iter1 = solid1.boundaries().iter(); + let shell0 = iter0.next().unwrap(); + let shell1 = iter1.next().unwrap(); + let [mut and_shell, _] = process_one_pair_of_shells(shell0, shell1, tol)?; + for shell in iter0 { + let [res, _] = process_one_pair_of_shells(&and_shell, shell, tol)?; + and_shell = res; + } + for shell in iter1 { + let [res, _] = process_one_pair_of_shells(&and_shell, shell, tol)?; + and_shell = res; + } + let boundaries = and_shell.connected_components(); + Some(Solid::new(boundaries)) +} + +/// OR operation between two solids. +pub fn or, S: ShapeOpsSurface>( + solid0: &Solid, + solid1: &Solid, + tol: f64, +) -> Option> { + let mut iter0 = solid0.boundaries().iter(); + let mut iter1 = solid1.boundaries().iter(); + let shell0 = iter0.next().unwrap(); + let shell1 = iter1.next().unwrap(); + let [_, mut or_shell] = process_one_pair_of_shells(shell0, shell1, tol)?; + for shell in iter0 { + let [_, res] = process_one_pair_of_shells(&or_shell, shell, tol)?; + or_shell = res; + } + for shell in iter1 { + let [_, res] = process_one_pair_of_shells(&or_shell, shell, tol)?; + or_shell = res; + } + let boundaries = or_shell.connected_components(); + Some(Solid::new(boundaries)) +} + +#[cfg(test)] +mod tests; diff --git a/vendor/truck-shapeops/src/transversal/integrate/tests.rs b/vendor/truck-shapeops/src/transversal/integrate/tests.rs new file mode 100644 index 00000000..b59d87ef --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/integrate/tests.rs @@ -0,0 +1,21 @@ +use truck_meshalgo::prelude::*; +use truck_modeling::*; + +#[test] +fn punched_cube() { + let v = builder::vertex(Point3::origin()); + let e = builder::tsweep(&v, Vector3::unit_x()); + let f = builder::tsweep(&e, Vector3::unit_y()); + let cube = builder::tsweep(&f, Vector3::unit_z()); + + let v = builder::vertex(Point3::new(0.5, 0.25, -0.5)); + let w = builder::rsweep(&v, Point3::new(0.5, 0.5, 0.0), Vector3::unit_z(), Rad(7.0)); + let f = builder::try_attach_plane(&[w]).unwrap(); + let mut cylinder = builder::tsweep(&f, Vector3::unit_z() * 2.0); + cylinder.not(); + let and = crate::and(&cube, &cylinder, 0.05).unwrap(); + + let poly = and.triangulation(0.01).to_polygon(); + let file = std::fs::File::create("punched-cube.obj").unwrap(); + obj::write(&poly, file).unwrap(); +} diff --git a/vendor/truck-shapeops/src/transversal/intersection_curve/mod.rs b/vendor/truck-shapeops/src/transversal/intersection_curve/mod.rs new file mode 100644 index 00000000..3988574d --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/intersection_curve/mod.rs @@ -0,0 +1,174 @@ +use truck_base::cgmath64::*; +use truck_geometry::prelude::*; +use truck_meshalgo::prelude::*; + +/// polyline base ntersection curve with parameter +#[derive(Debug, Clone, derive_more::Deref, derive_more::DerefMut)] +pub struct IntersectionCurveWithParameters { + #[deref] + #[deref_mut] + ic: IntersectionCurve, S>, + params0: PolylineCurve, + params1: PolylineCurve, +} + +impl From> for IntersectionCurve, S> { + fn from(a: IntersectionCurveWithParameters) -> Self { a.ic } +} + +impl IntersectionCurveWithParameters +where S: ParametricSurface3D + SearchNearestParameter +{ + pub fn try_new( + surface0: S, + surface1: S, + poly: PolylineCurve, + tol: f64, + ) -> Option { + let mut polyline = PolylineCurve(Vec::new()); + let mut params0 = PolylineCurve(Vec::new()); + let mut params1 = PolylineCurve(Vec::new()); + for p in poly.windows(2) { + let n = (p[1] - p[0]).normalize(); + let (q, p0, p1) = double_projection(&surface0, None, &surface1, None, p[0], n, 100)?; + polyline.push(q); + params0.push(p0); + params1.push(p1); + } + let (q, p0, p1) = if poly[0].near(&poly[poly.len() - 1]) { + (polyline[0], params0[0], params1[0]) + } else { + let n = (poly[poly.len() - 1] - poly[poly.len() - 2]).normalize(); + double_projection( + &surface0, + None, + &surface1, + None, + poly[poly.len() - 1], + n, + 100, + )? + }; + polyline.push(q); + params0.push(p0); + params1.push(p1); + Some(Self { + ic: IntersectionCurve::new_unchecked( + Box::new(surface0), + Box::new(surface1), + polyline, + tol, + ), + params0, + params1, + }) + } +} + +impl ParametricCurve for IntersectionCurveWithParameters +where S: ParametricSurface3D + SearchNearestParameter +{ + type Point = Point3; + type Vector = Vector3; + #[inline(always)] + fn subs(&self, t: f64) -> Point3 { self.ic.subs(t) } + #[inline(always)] + fn der(&self, t: f64) -> Vector3 { self.ic.der(t) } + #[inline(always)] + fn der2(&self, t: f64) -> Vector3 { self.ic.der2(t) } + #[inline(always)] + fn parameter_range(&self) -> ParameterRange { self.ic.parameter_range() } +} + +impl BoundedCurve for IntersectionCurveWithParameters where S: ParametricSurface3D + SearchNearestParameter {} + +impl ParameterDivision1D for IntersectionCurveWithParameters +where S: ParametricSurface3D + SearchNearestParameter +{ + type Point = Point3; + #[inline(always)] + fn parameter_division(&self, range: (f64, f64), tol: f64) -> (Vec, Vec) { + self.ic.parameter_division(range, tol) + } +} + +impl Cut for IntersectionCurveWithParameters +where S: ParametricSurface3D + SearchNearestParameter +{ + #[inline(always)] + fn cut(&mut self, t: f64) -> Self { + Self { + ic: self.ic.cut(t), + params0: self.params0.cut(t), + params1: self.params1.cut(t), + } + } +} + +impl Invertible for IntersectionCurveWithParameters { + fn invert(&mut self) { + self.ic.invert(); + self.params0.invert(); + self.params1.invert(); + } +} + +impl SearchParameter for IntersectionCurveWithParameters +where S: ParametricSurface3D + SearchNearestParameter +{ + type Point = Point3; + fn search_parameter>( + &self, + point: Point3, + hint: H, + trials: usize, + ) -> Option { + self.ic.search_parameter(point, hint, trials) + } +} + +impl SearchNearestParameter for IntersectionCurveWithParameters +where S: ParametricSurface3D + SearchNearestParameter +{ + type Point = Point3; + fn search_nearest_parameter>( + &self, + point: Point3, + hint: H, + trials: usize, + ) -> Option { + self.ic.search_nearest_parameter(point, hint, trials) + } +} + +pub fn intersection_curves( + surface0: S, + polygon0: &PolygonMesh, + surface1: S, + polygon1: &PolygonMesh, + tol: f64, +) -> Vec<( + PolylineCurve, + Option>, +)> +where + S: ParametricSurface3D + SearchNearestParameter, +{ + let interferences = polygon0.extract_interference(polygon1); + let polylines = super::polyline_construction::construct_polylines(&interferences); + polylines + .into_iter() + .map(|polyline| { + let curve = IntersectionCurveWithParameters::try_new( + surface0.clone(), + surface1.clone(), + polyline.clone(), + tol, + ); + (polyline, curve) + }) + .collect() +} + +#[cfg(test)] +mod tests; diff --git a/vendor/truck-shapeops/src/transversal/intersection_curve/tests.rs b/vendor/truck-shapeops/src/transversal/intersection_curve/tests.rs new file mode 100644 index 00000000..1368b773 --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/intersection_curve/tests.rs @@ -0,0 +1,173 @@ +use super::*; +use std::f64::consts::PI; + +mod double_projection_tests { + use super::*; + + fn create_axis(n: Vector3) -> (Vector3, Vector3) { + let idx = if n[0].abs() < n[1].abs() { 0 } else { 1 }; + let idx = if n[idx].abs() < n[2].abs() { idx } else { 2 }; + let mut e = Vector3::zero(); + e[idx] = 1.0; + let x = n.cross(e).normalize(); + (x, n.cross(x)) + } + + fn exec_plane_case() { + let c0 = Point3::new( + 2.0 * rand::random::() - 1.0, + 2.0 * rand::random::() - 1.0, + 2.0 * rand::random::() - 1.0, + ); + let n0 = Vector3::new( + 2.0 * rand::random::() - 1.0, + 2.0 * rand::random::() - 1.0, + 2.0 * rand::random::() - 1.0, + ) + .normalize(); + let (x, y) = create_axis(n0); + let plane0 = Plane::new(c0, c0 + x, c0 + y); + let c1 = Point3::new( + 2.0 * rand::random::() - 1.0, + 2.0 * rand::random::() - 1.0, + 2.0 * rand::random::() - 1.0, + ); + let n1 = Vector3::new( + 2.0 * rand::random::() - 1.0, + 2.0 * rand::random::() - 1.0, + 2.0 * rand::random::() - 1.0, + ) + .normalize(); + let (x, y) = create_axis(n1); + let plane1 = Plane::new(c1, c1 + x, c1 + y); + let n = n0.cross(n1).normalize(); + let mut o = None; + for i in 0..10 { + let t = i as f64; + let p = Point3::origin() + t * n; + let (q, p0, p1) = double_projection(&plane0, None, &plane1, None, p, n, 100) + .unwrap_or_else(|| panic!("plane0: {plane0:?}\nplane1: {plane1:?}\n p: {p:?}")); + assert_near!(q, plane0.subs(p0.x, p0.y)); + assert_near!(q, plane1.subs(p1.x, p1.y)); + if let Some(o) = o { + assert_near!(q.distance2(o), t * t); + } else { + o = Some(q); + } + } + } + + #[test] + fn plane_case() { (0..100).for_each(|_| exec_plane_case()); } + + #[test] + fn sphere_case() { + let sphere0 = Sphere::new(Point3::new(0.0, 0.0, 1.0), f64::sqrt(2.0)); + let sphere1 = Sphere::new(Point3::new(0.0, 0.0, -1.0), f64::sqrt(2.0)); + for _ in 0..100 { + let t = 2.0 * PI * rand::random::(); + let r = 2.0 * rand::random::(); + let p = Point3::new(r * f64::cos(t), r * f64::sin(t), 0.0); + let n = Vector3::new(-f64::sin(t), f64::cos(t), 0.0); + let (q, p0, p1) = double_projection(&sphere0, None, &sphere1, None, p, n, 100) + .unwrap_or_else(|| panic!("p: {p:?}")); + assert_near!(q, sphere0.subs(p0.x, p0.y)); + assert_near!(q, sphere1.subs(p1.x, p1.y)); + assert_near!(q, Point3::new(f64::cos(t), f64::sin(t), 0.0)); + } + } +} + +#[test] +fn intersection_curve_sphere_case() { + let sphere0 = Sphere::new(Point3::new(0.0, 0.0, 1.0), f64::sqrt(2.0)); + let sphere1 = Sphere::new(Point3::new(0.0, 0.0, -1.0), f64::sqrt(2.0)); + const M: usize = 5; + let polyline = (0..=M) + .map(|i| { + let t = 2.0 * PI * i as f64 / M as f64; + Point3::new(0.8 * f64::cos(t), 0.8 * f64::sin(t), 0.0) + }) + .collect::>(); + let curve: IntersectionCurve<_, _> = + IntersectionCurveWithParameters::try_new(sphere0, sphere1, polyline, 0.5) + .unwrap() + .into(); + + const N: usize = 100; + let mut sum = 0.0; + let (t0, t1) = curve.range_tuple(); + for i in 0..N { + let t = t0 + (t1 - t0) * i as f64 / N as f64; + let pt = curve.subs(t); + assert_near!(pt.distance(Point3::origin()), 1.0); + let vec = curve.der(t); + assert!(pt.dot(vec).so_small(), "{i} {t} {vec:?}"); + assert!(vec[2].so_small()); + sum += vec.magnitude() * (t1 - t0) / N as f64; + } + assert!(f64::abs(sum - 2.0 * PI) < 0.1, "{}", sum); + + let theta = 2.0 * PI * rand::random::(); + let pt = Point3::new(f64::cos(theta), f64::sin(theta), 0.0); + let t = curve.search_parameter(pt, None, 1).unwrap(); + assert!(curve.subs(t).near(&pt)); + let pt = Point3::new(1.1 * f64::cos(theta), 1.1 * f64::sin(theta), 0.0); + assert!(curve.search_parameter(pt, None, 1).is_none()); + + let mut curve0 = curve.clone(); + let curve1 = curve0.cut(2.5); + assert_near!(curve0.front(), curve.front()); + assert_near!(curve0.back(), curve.subs(2.5)); + assert_near!(curve1.front(), curve.subs(2.5)); + assert_near!(curve1.back(), curve.back()); + let mut curve0 = curve.clone(); + let curve1 = curve0.cut(2.0); + assert_near!(curve0.front(), curve.front()); + assert_near!(curve0.back(), curve.subs(2.0)); + assert_near!(curve1.front(), curve.subs(2.0)); + assert_near!(curve1.back(), curve.back()); +} + +#[test] +fn collide_parabola() { + const TOL: f64 = 0.05; + + // define surfaces + #[rustfmt::skip] + let ctrl0 = vec![ + vec![Point3::new(-1.0, -1.0, 3.0), Point3::new(-1.0, 0.0, -1.0), Point3::new(-1.0, 1.0, 3.0)], + vec![Point3::new(0.0, -1.0, -1.0), Point3::new(0.0, 0.0, -5.0), Point3::new(0.0, 1.0, -1.0)], + vec![Point3::new(1.0, -1.0, 3.0), Point3::new(1.0, 0.0, -1.0), Point3::new(1.0, 1.0, 3.0)], + ]; + #[rustfmt::skip] + let ctrl1 = vec![ + vec![Point3::new(-1.0, -1.0, -3.0), Point3::new(-1.0, 0.0, 1.0), Point3::new(-1.0, 1.0, -3.0)], + vec![Point3::new(0.0, -1.0, 1.0), Point3::new(0.0, 0.0, 5.0), Point3::new(0.0, 1.0, 1.0)], + vec![Point3::new(1.0, -1.0, -3.0), Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 1.0, -3.0)], + ]; + let surface0 = BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl0); + let surface1 = BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl1); + + // meshing surface + let instant = std::time::Instant::now(); + let polygon0 = StructuredMesh::from_surface(&surface0, surface0.range_tuple(), TOL).destruct(); + let polygon1 = StructuredMesh::from_surface(&surface1, surface1.range_tuple(), TOL).destruct(); + println!("Meshing Surfaces: {}s", instant.elapsed().as_secs_f64()); + // extract intersection curves + let instant = std::time::Instant::now(); + let curves = intersection_curves(surface0, &polygon0, surface1, &polygon1, TOL); + println!( + "Extracting Intersection: {}s", + instant.elapsed().as_secs_f64() + ); + assert_eq!(curves.len(), 1); + let curve = curves[0].1.clone().unwrap(); + const N: usize = 100; + for i in 0..N { + let t1 = curve.range_tuple().1; + let t = t1 * i as f64 / N as f64; + let pt = curve.subs(t); + assert_near!(pt.distance(Point3::origin()) * 0.5, f64::sqrt(0.5) * 0.5); + } +} diff --git a/vendor/truck-shapeops/src/transversal/loops_store/mod.rs b/vendor/truck-shapeops/src/transversal/loops_store/mod.rs new file mode 100644 index 00000000..7958997e --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/loops_store/mod.rs @@ -0,0 +1,579 @@ +#![allow(clippy::many_single_char_names)] + +use super::*; +use rustc_hash::FxHashMap as HashMap; +use truck_base::cgmath64::*; +use truck_geometry::prelude::IntersectionCurve; +use truck_meshalgo::prelude::*; +use truck_topology::{Vertex, *}; + +type PolylineCurve = truck_meshalgo::prelude::PolylineCurve; + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum ShapesOpStatus { + Unknown, + And, + Or, +} + +impl ShapesOpStatus { + fn not(self) -> Self { + match self { + Self::Unknown => Self::Unknown, + Self::And => Self::Or, + Self::Or => Self::And, + } + } +} + +#[derive(Clone, Debug)] +pub struct BoundaryWire { + wire: Wire, + status: ShapesOpStatus, +} + +impl BoundaryWire { + #[inline(always)] + pub fn new(wire: Wire, status: ShapesOpStatus) -> Self { Self { wire, status } } + #[inline(always)] + pub fn status(&self) -> ShapesOpStatus { self.status } + #[inline(always)] + pub fn invert(&mut self) { + self.wire.invert(); + self.status = self.status.not(); + } + #[inline(always)] + pub fn inverse(&self) -> Self { + Self { + wire: self.wire.inverse(), + status: self.status.not(), + } + } +} + +impl ShapesOpStatus { + fn from_is_curve(curve: &IntersectionCurve) -> Option + where + C: ParametricCurve3D + BoundedCurve, + S: ParametricSurface3D + SearchNearestParameter, { + let (t0, t1) = curve.range_tuple(); + let t = (t0 + t1) / 2.0; + let (_, pt0, pt1) = curve.search_triple(t)?; + let der = curve.leader().der(t); + let normal0 = curve.surface0().normal(pt0[0], pt0[1]); + let normal1 = curve.surface1().normal(pt1[0], pt1[1]); + match normal0.cross(der).dot(normal1) > 0.0 { + true => Some(ShapesOpStatus::Or), + false => Some(ShapesOpStatus::And), + } + } +} + +impl std::ops::Deref for BoundaryWire { + type Target = Wire; + #[inline(always)] + fn deref(&self) -> &Self::Target { &self.wire } +} + +impl std::ops::DerefMut for BoundaryWire { + #[inline(always)] + fn deref_mut(&mut self) -> &mut Self::Target { &mut self.wire } +} + +#[derive(Clone, Debug)] +pub struct Loops(Vec>); +#[derive(Clone, Debug)] +pub struct LoopsStore(Vec>); + +impl std::ops::Deref for Loops { + type Target = Vec>; + #[inline(always)] + fn deref(&self) -> &Self::Target { &self.0 } +} + +impl std::ops::DerefMut for Loops { + #[inline(always)] + fn deref_mut(&mut self) -> &mut Self::Target { &mut self.0 } +} + +impl std::ops::Deref for LoopsStore { + type Target = Vec>; + #[inline(always)] + fn deref(&self) -> &Self::Target { &self.0 } +} + +impl std::ops::DerefMut for LoopsStore { + #[inline(always)] + fn deref_mut(&mut self) -> &mut Self::Target { &mut self.0 } +} + +impl FromIterator> for Loops { + #[inline(always)] + fn from_iter>>(iter: I) -> Self { + Self(Vec::from_iter(iter)) + } +} + +impl<'a, P, C, S> From<&'a Face> for Loops { + #[inline(always)] + fn from(face: &'a Face) -> Loops { + face.absolute_boundaries() + .iter() + .map(|wire| BoundaryWire::new(wire.clone(), ShapesOpStatus::Unknown)) + .collect() + } +} + +impl<'a, P: 'a, C: 'a, S: 'a> FromIterator<&'a Face> for LoopsStore { + fn from_iter>>(iter: I) -> Self { + Self(iter.into_iter().map(Loops::from).collect()) + } +} + +impl<'a, P, C> IntoIterator for &'a LoopsStore { + type Item = <&'a Vec> as IntoIterator>::Item; + type IntoIter = <&'a Vec> as IntoIterator>::IntoIter; + fn into_iter(self) -> Self::IntoIter { self.0.iter() } +} + +#[derive(Clone, Debug, Copy, PartialEq)] +enum ParameterKind { + Front, + Back, + Inner(f64), +} + +impl ParameterKind { + fn try_new(t: f64, (t0, t1): (f64, f64)) -> Option { + if t0.near(&t) { + Some(ParameterKind::Front) + } else if t1.near(&t) { + Some(ParameterKind::Back) + } else if t0 < t && t < t1 { + Some(ParameterKind::Inner(t)) + } else { + None + } + } +} + +impl Loops { + fn search_parameter(&self, pt: P) -> Option<(usize, usize, ParameterKind)> + where C: BoundedCurve + SearchParameter { + self.iter() + .enumerate() + .flat_map(move |(i, wire)| wire.iter().enumerate().map(move |(j, edge)| (i, j, edge))) + .find_map(|(i, j, edge)| { + let curve = edge.curve(); + curve.search_parameter(pt, None, 1).and_then(|t| { + let kind = ParameterKind::try_new(t, curve.range_tuple())?; + Some((i, j, kind)) + }) + }) + } + + fn change_vertex( + &mut self, + old_vertex: &Vertex

, + new_vertex: &Vertex

, + emap: &mut HashMap, Edge>, + ) { + self.iter_mut() + .flat_map(|wire| wire.iter_mut()) + .for_each(|edge| { + let mut new_edge = if edge.absolute_front() == old_vertex { + emap.entry(edge.id()).or_insert_with(|| { + Edge::new(new_vertex, edge.absolute_back(), edge.curve()) + }) + } else if edge.absolute_back() == old_vertex { + emap.entry(edge.id()).or_insert_with(|| { + Edge::new(edge.absolute_front(), new_vertex, edge.curve()) + }) + } else { + return; + } + .clone(); + if !edge.orientation() { + new_edge.invert(); + } + // Remove the edge from the HashMap when it is no longer there because ID reassignment will occur. + if edge.count() == 1 { + emap.remove(&edge.id()); + } + *edge = new_edge; + }) + } + + fn swap_edge_into_wire(&mut self, edge_id: EdgeID, new_wire: &Wire) { + self.iter_mut().for_each(|wire| { + let mut iter = wire.iter().enumerate(); + if let Some((idx, edge)) = iter.find(|(_, edge)| edge.id() == edge_id) { + if edge.orientation() { + wire.swap_edge_into_wire(idx, new_wire.clone()); + } else { + wire.swap_edge_into_wire(idx, new_wire.inverse()); + } + } + }); + } + + #[inline(always)] + fn add_independent_loop(&mut self, r#loop: BoundaryWire) { + self.push(r#loop.inverse()); + self.push(r#loop); + } + + fn add_edge( + &mut self, + edge0: Edge, + status: ShapesOpStatus, + ) -> [Option<(usize, usize)>; 2] { + let a = self.iter().enumerate().find_map(|(i, wire)| { + wire.iter().enumerate().find_map(|(j, edge)| { + if edge.front() == edge0.back() { + Some((i, j)) + } else { + None + } + }) + }); + let b = self.iter().enumerate().find_map(|(i, wire)| { + wire.iter().enumerate().find_map(|(j, edge)| { + if edge.front() == edge0.front() { + Some((i, j)) + } else { + None + } + }) + }); + if let Some((wire_index0, edge_index0)) = a { + self[wire_index0].rotate_left(edge_index0); + self[wire_index0].push_front(edge0.clone()); + self[wire_index0].push_back(edge0.inverse()); + } + match (a, b) { + (Some((wire_index0, edge_index0)), Some((wire_index1, edge_index1))) => { + if wire_index0 == wire_index1 { + let len = self[wire_index0].len() - 2; + let edge_index1 = (len + edge_index1 - edge_index0) % len + 1; + let new_wire = self[wire_index0].split_off(edge_index1); + self[wire_index0].status = status; + self.push(BoundaryWire::new(new_wire, status.not())); + } else { + let mut new_wire0 = self[wire_index1].clone(); + let mut new_wire1 = new_wire0.split_off(edge_index1); + new_wire0.append(&mut self[wire_index0]); + new_wire0.append(&mut new_wire1); + self[wire_index0] = new_wire0; + self.swap_remove(wire_index1); + } + } + (None, Some((wire_index1, edge_index1))) => { + self[wire_index1].rotate_left(edge_index1); + self[wire_index1].push_front(edge0.inverse()); + self[wire_index1].push_back(edge0); + } + (None, None) => self.push(BoundaryWire::new( + vec![edge0.inverse(), edge0].into(), + ShapesOpStatus::Unknown, + )), + _ => {} + } + [a, b] + } +} + +impl LoopsStore { + #[inline(always)] + fn change_vertex( + &mut self, + old_vertex: &Vertex

, + new_vertex: &Vertex

, + emap: &mut HashMap, Edge>, + ) { + self.iter_mut() + .for_each(|loops| loops.change_vertex(old_vertex, new_vertex, emap)); + } + + #[inline(always)] + fn swap_edge_into_wire(&mut self, edge_id: EdgeID, new_wire: &Wire) { + self.iter_mut() + .for_each(|loops| loops.swap_edge_into_wire(edge_id, new_wire)) + } + + fn add_polygon_vertex( + &mut self, + loops_index: usize, + v: &Vertex

, + emap: &mut HashMap, Edge>, + ) -> Option<(usize, usize, ParameterKind)> + where + C: Cut + SearchParameter, + { + let pt = v.point(); + let (wire_index, edge_index, kind) = self[loops_index].search_parameter(pt)?; + match kind { + ParameterKind::Front => { + let old_vertex = self[loops_index][wire_index][edge_index] + .absolute_front() + .clone(); + self.change_vertex(&old_vertex, v, emap); + } + ParameterKind::Back => { + let old_vertex = self[loops_index][wire_index][edge_index] + .absolute_back() + .clone(); + self.change_vertex(&old_vertex, v, emap); + } + ParameterKind::Inner(t) => { + let edge = self[loops_index][wire_index][edge_index].absolute_clone(); + let edge_id = edge.id(); + let (edge0, edge1) = edge.cut_with_parameter(v, t)?; + let new_wire: Wire<_, _> = vec![edge0, edge1].into(); + self.swap_edge_into_wire(edge_id, &new_wire); + } + } + Some((wire_index, edge_index, kind)) + } +} + +impl LoopsStore { + fn add_geom_vertex( + &mut self, + (loops_index, wire_index, edge_index): (usize, usize, usize), + v: &Vertex, + kind: ParameterKind, + another_surface: &S, + emap: &mut HashMap, Edge>, + ) -> Option<()> + where + C: Cut + SearchNearestParameter, + S: ParametricSurface3D + SearchNearestParameter, + { + match kind { + ParameterKind::Front => { + let old_vertex = self[loops_index][wire_index][edge_index] + .absolute_front() + .clone(); + v.set_point(old_vertex.point()); + self.change_vertex(&old_vertex, v, emap); + } + ParameterKind::Back => { + let old_vertex = self[loops_index][wire_index][edge_index] + .absolute_back() + .clone(); + v.set_point(old_vertex.point()); + self.change_vertex(&old_vertex, v, emap); + } + ParameterKind::Inner(_) => { + let curve = self[loops_index][wire_index][edge_index].curve(); + let (pt, t, _) = + curve_surface_projection(&curve, None, another_surface, None, v.point(), 100)?; + v.set_point(pt); + let edge = self[loops_index][wire_index][edge_index].absolute_clone(); + let edge_id = edge.id(); + let (edge0, edge1) = edge.cut_with_parameter(v, t)?; + let new_wire: Wire<_, _> = vec![edge0, edge1].into(); + self.swap_edge_into_wire(edge_id, &new_wire); + } + } + Some(()) + } +} + +fn curve_surface_projection( + curve: &C, + curve_hint: Option, + surface: &S, + surface_hint: Option<(f64, f64)>, + point: Point3, + trials: usize, +) -> Option<(Point3, f64, Point2)> +where + C: ParametricCurve3D + SearchNearestParameter, + S: ParametricSurface3D + SearchNearestParameter, +{ + if trials == 0 { + return None; + } + let t = curve.search_nearest_parameter(point, curve_hint, 10)?; + let pt0 = curve.subs(t); + let (u, v) = surface.search_nearest_parameter(point, surface_hint, 10)?; + let pt1 = surface.subs(u, v); + if point.near(&pt0) && point.near(&pt1) && pt0.near(&pt1) { + Some((point, t, Point2::new(u, v))) + } else { + let l = curve.der(t); + let n = surface.normal(u, v); + let t0 = (pt1 - pt0).dot(n) / l.dot(n); + curve_surface_projection( + curve, + Some(t), + surface, + Some((u, v)), + pt0 + t0 * l, + trials - 1, + ) + } +} + +fn create_independent_loop(mut poly_curve0: C) -> Wire +where + C: Cut, + D: From, { + let (t0, t1) = poly_curve0.range_tuple(); + let t = (t0 + t1) / 2.0; + let poly_curve1 = poly_curve0.cut(t); + let v0 = Vertex::new(poly_curve0.front()); + let v1 = Vertex::new(poly_curve1.front()); + let edge0 = Edge::new(&v0, &v1, poly_curve0.into()); + let edge1 = Edge::new(&v1, &v0, poly_curve1.into()); + vec![edge0, edge1].into() +} + +#[allow(dead_code)] +pub struct LoopsStoreQuadruple { + pub geom_loops_store0: LoopsStore, + pub poly_loops_store0: LoopsStore, + pub geom_loops_store1: LoopsStore, + pub poly_loops_store1: LoopsStore, +} + +pub fn create_loops_stores( + geom_shell0: &Shell, + poly_shell0: &Shell>, + geom_shell1: &Shell, + poly_shell1: &Shell>, + tol: f64, +) -> Option> +where + C: SearchNearestParameter + + SearchParameter + + Cut + + From>, + S: ParametricSurface3D + SearchNearestParameter, +{ + let mut geom_loops_store0: LoopsStore<_, _> = geom_shell0.face_iter().collect(); + let mut poly_loops_store0: LoopsStore<_, _> = poly_shell0.face_iter().collect(); + let mut geom_loops_store1: LoopsStore<_, _> = geom_shell1.face_iter().collect(); + let mut poly_loops_store1: LoopsStore<_, _> = poly_shell1.face_iter().collect(); + let store0_len = geom_loops_store0.len(); + let store1_len = geom_loops_store1.len(); + (0..store0_len) + .flat_map(move |i| (0..store1_len).map(move |j| (i, j))) + .try_for_each(|(face_index0, face_index1)| { + let ori0 = geom_shell0[face_index0].orientation(); + let ori1 = geom_shell1[face_index1].orientation(); + let surface0 = geom_shell0[face_index0].surface(); + let surface1 = geom_shell1[face_index1].surface(); + let polygon0 = crate::tr!(poly_shell0[face_index0].surface(), "polygon0"); + let polygon1 = crate::tr!(poly_shell1[face_index1].surface(), "polygon1"); + intersection_curve::intersection_curves( + surface0.clone(), + &polygon0, + surface1.clone(), + &polygon1, + tol, + ) + .into_iter() + .try_for_each(|(polyline, intersection_curve)| { + let mut intersection_curve = crate::tr!(intersection_curve, "intersection curve").into(); + let status = crate::tr!(ShapesOpStatus::from_is_curve(&intersection_curve), "status"); + let (status0, status1) = match (ori0, ori1) { + (true, true) => (status, status.not()), + (true, false) => (status.not(), status.not()), + (false, true) => (status, status), + (false, false) => (status.not(), status), + }; + if polyline.front().near(&polyline.back()) { + let poly_wire = create_independent_loop(polyline); + poly_loops_store0[face_index0] + .add_independent_loop(BoundaryWire::new(poly_wire.clone(), status0)); + poly_loops_store1[face_index1] + .add_independent_loop(BoundaryWire::new(poly_wire, status1)); + let geom_wire = create_independent_loop(intersection_curve); + geom_loops_store0[face_index0] + .add_independent_loop(BoundaryWire::new(geom_wire.clone(), status0)); + geom_loops_store1[face_index1] + .add_independent_loop(BoundaryWire::new(geom_wire, status1)); + } else { + let pv0 = Vertex::new(polyline.front()); + let pv1 = Vertex::new(polyline.back()); + let gv0 = Vertex::new(polyline.front()); + let gv1 = Vertex::new(polyline.back()); + let mut pemap0 = HashMap::default(); + let mut pemap1 = HashMap::default(); + let mut gemap0 = HashMap::default(); + let mut gemap1 = HashMap::default(); + let idx00 = + poly_loops_store0.add_polygon_vertex(face_index0, &pv0, &mut pemap0); + if let Some((wire_index, edge_index, kind)) = idx00 { + crate::tr!(geom_loops_store0.add_geom_vertex( + (face_index0, wire_index, edge_index), + &gv0, + kind, + &surface1, + &mut gemap0, + ), "add_geom_vertex 1"); + let polyline = intersection_curve.editor().leader; + *polyline.first_mut().unwrap() = gv0.point(); + } + let idx01 = + poly_loops_store0.add_polygon_vertex(face_index0, &pv1, &mut pemap1); + if let Some((wire_index, edge_index, kind)) = idx01 { + crate::tr!(geom_loops_store0.add_geom_vertex( + (face_index0, wire_index, edge_index), + &gv1, + kind, + &surface1, + &mut gemap1, + ), "add_geom_vertex 2"); + let polyline = intersection_curve.editor().leader; + *polyline.last_mut().unwrap() = gv1.point(); + } + let idx10 = + poly_loops_store1.add_polygon_vertex(face_index1, &pv0, &mut pemap0); + if let Some((wire_index, edge_index, kind)) = idx10 { + crate::tr!(geom_loops_store1.add_geom_vertex( + (face_index1, wire_index, edge_index), + &gv0, + kind, + &surface0, + &mut gemap0, + ), "add_geom_vertex 3"); + let polyline = intersection_curve.editor().leader; + *polyline.first_mut().unwrap() = gv0.point(); + } + let idx11 = + poly_loops_store1.add_polygon_vertex(face_index1, &pv1, &mut pemap1); + if let Some((wire_index, edge_index, kind)) = idx11 { + crate::tr!(geom_loops_store1.add_geom_vertex( + (face_index1, wire_index, edge_index), + &gv1, + kind, + &surface0, + &mut gemap1, + ), "add_geom_vertex 4"); + let polyline = intersection_curve.editor().leader; + *polyline.last_mut().unwrap() = gv1.point(); + } + let pedge = Edge::new(&pv0, &pv1, polyline); + let gedge = Edge::new(&gv0, &gv1, intersection_curve.into()); + poly_loops_store0[face_index0].add_edge(pedge.clone(), status0); + geom_loops_store0[face_index0].add_edge(gedge.clone(), status0); + poly_loops_store1[face_index1].add_edge(pedge, status1); + geom_loops_store1[face_index1].add_edge(gedge, status1); + } + Some(()) + }) + })?; + Some(LoopsStoreQuadruple { + geom_loops_store0, + poly_loops_store0, + geom_loops_store1, + poly_loops_store1, + }) +} + +#[cfg(test)] +mod tests; diff --git a/vendor/truck-shapeops/src/transversal/loops_store/tests.rs b/vendor/truck-shapeops/src/transversal/loops_store/tests.rs new file mode 100644 index 00000000..defbad19 --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/loops_store/tests.rs @@ -0,0 +1,776 @@ +use super::*; +use std::fmt::Debug; +use truck_geometry::prelude::*; +const TOL: f64 = 0.01; + +crate::impl_from!( + NurbsCurve, + IntersectionCurve> +); +type AlternativeIntersection = crate::alternative::Alternative< + NurbsCurve, + IntersectionCurve>, +>; + +struct DebugDisplay<'a, T, Format> { + entity: &'a T, + format: Format, +} + +impl<'a, P: Debug, C: Debug> Debug for DebugDisplay<'a, Loops, WireDisplayFormat> { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_tuple("Loops") + .field( + &self + .entity + .0 + .iter() + .map(|wire| wire.display(self.format)) + .collect::>(), + ) + .finish() + } +} + +impl<'a, P: Debug, C: Debug> Debug for DebugDisplay<'a, LoopsStore, WireDisplayFormat> { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_list() + .entries(self.entity.0.iter().map(|loops| DebugDisplay { + entity: loops, + format: self.format, + })) + .finish() + } +} + +impl LoopsStore { + fn display(&self, format: WireDisplayFormat) -> DebugDisplay<'_, Self, WireDisplayFormat> { + DebugDisplay { + entity: self, + format, + } + } +} + +fn parabola_surfaces() -> (BSplineSurface, BSplineSurface) { + // define surfaces + #[rustfmt::skip] + let ctrl0 = vec![ + vec![Point3::new(-1.0, -1.0, 3.0), Point3::new(-1.0, 0.0, -1.0), Point3::new(-1.0, 1.0, 3.0)], + vec![Point3::new(0.0, -1.0, -1.0), Point3::new(0.0, 0.0, -5.0), Point3::new(0.0, 1.0, -1.0)], + vec![Point3::new(1.0, -1.0, 3.0), Point3::new(1.0, 0.0, -1.0), Point3::new(1.0, 1.0, 3.0)], + ]; + #[rustfmt::skip] + let ctrl1 = vec![ + vec![Point3::new(-1.0, -1.0, -3.0), Point3::new(-1.0, 0.0, 1.0), Point3::new(-1.0, 1.0, -3.0)], + vec![Point3::new(0.0, -1.0, 1.0), Point3::new(0.0, 0.0, 5.0), Point3::new(0.0, 1.0, 1.0)], + vec![Point3::new(1.0, -1.0, -3.0), Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 1.0, -3.0)], + ]; + ( + BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl0), + BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl1), + ) +} + +#[test] +fn independent_intersection() { + // prepare geoetries + let arc00: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, 1.0, 0.0, 0.0), + Vector4::new(-1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc01 = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(-1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, -1.0, 0.0, 0.0), + Vector4::new(1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc10: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, -1.0, 1.0), + Vector4::new(0.0, 1.0, 0.0, 0.0), + Vector4::new(-1.0, 0.0, -1.0, 1.0), + ], + )) + .into(); + let arc11 = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(-1.0, 0.0, -1.0, 1.0), + Vector4::new(0.0, -1.0, 0.0, 0.0), + Vector4::new(1.0, 0.0, -1.0, 1.0), + ], + )) + .into(); + let (surface0, surface1) = parabola_surfaces(); + + // prepare topologies + let v00 = Vertex::new(Point3::new(1.0, 0.0, 1.0)); + let v01 = Vertex::new(Point3::new(-1.0, 0.0, 1.0)); + let v10 = Vertex::new(Point3::new(1.0, 0.0, -1.0)); + let v11 = Vertex::new(Point3::new(-1.0, 0.0, -1.0)); + let wire0: Wire<_, _> = vec![Edge::new(&v00, &v01, arc00), Edge::new(&v01, &v00, arc01)].into(); + let wire1: Wire<_, _> = vec![Edge::new(&v10, &v11, arc10), Edge::new(&v11, &v10, arc11)].into(); + let geom_shell0: Shell<_, _, _> = vec![Face::new(vec![wire0], surface0).inverse()].into(); + let geom_shell1: Shell<_, _, _> = vec![Face::new(vec![wire1], surface1)].into(); + let poly_shell0 = geom_shell0.triangulation(TOL); + let poly_shell1 = geom_shell1.triangulation(TOL); + + // exec create loops stores! + let LoopsStoreQuadruple { + geom_loops_store0, + geom_loops_store1, + .. + } = create_loops_stores(&geom_shell0, &poly_shell0, &geom_shell1, &poly_shell1, TOL).unwrap(); + + // check the topology + let vertex_format = VertexDisplayFormat::AsPoint; + let edge_id_format = EdgeDisplayFormat::VerticesTupleAndID { vertex_format }; + let wire_id_format = WireDisplayFormat::EdgesListTuple { + edge_format: edge_id_format, + }; + let edge_geom_format = EdgeDisplayFormat::VerticesTupleAndCurve { vertex_format }; + let wire_geom_format = WireDisplayFormat::EdgesListTuple { + edge_format: edge_geom_format, + }; + assert_eq!(geom_loops_store0.len(), 1); + assert_eq!(geom_loops_store0[0].len(), 3); + assert_eq!(geom_loops_store0[0][0].len(), 2); + assert_eq!(geom_loops_store0[0][1].len(), 2); + assert_eq!(geom_loops_store0[0][2].len(), 2); + assert!( + geom_loops_store0[0][0].is_closed(), + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert!( + geom_loops_store0[0][1].is_closed(), + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert!( + geom_loops_store0[0][2].is_closed(), + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert!( + geom_loops_store0[0][0].is_geometric_consistent(), + "{:?}", + geom_loops_store0.display(wire_geom_format) + ); + assert!( + geom_loops_store0[0][1].is_geometric_consistent(), + "{:?}", + geom_loops_store0.display(wire_geom_format) + ); + assert!( + geom_loops_store0[0][2].is_geometric_consistent(), + "{:?}", + geom_loops_store0.display(wire_geom_format) + ); + assert_eq!(geom_loops_store1.len(), 1); + assert_eq!(geom_loops_store1[0].len(), 3); + assert_eq!(geom_loops_store1[0][0].len(), 2); + assert_eq!(geom_loops_store1[0][1].len(), 2); + assert_eq!(geom_loops_store1[0][2].len(), 2); + assert!( + geom_loops_store1[0][0].is_closed(), + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert!( + geom_loops_store1[0][1].is_closed(), + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert!( + geom_loops_store1[0][2].is_closed(), + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert!( + geom_loops_store1[0][0].is_geometric_consistent(), + "{:?}", + geom_loops_store1.display(wire_geom_format) + ); + assert!( + geom_loops_store1[0][1].is_geometric_consistent(), + "{:?}", + geom_loops_store1.display(wire_geom_format) + ); + assert!( + geom_loops_store1[0][2].is_geometric_consistent(), + "{:?}", + geom_loops_store1.display(wire_geom_format) + ); + + // check the boundary status + let mut flags = [true; 3]; + for i in 0..3 { + let bw = &geom_loops_store0[0][i]; + match bw.status() { + ShapesOpStatus::Unknown => { + let curve = bw[0].oriented_curve(); + assert_near!(curve.subs(0.0)[2], 1.0); + assert!(flags[0]); + flags[0] = false; + } + ShapesOpStatus::And => { + let curve = bw[0].oriented_curve(); + let pt = curve.subs(0.5) - Point3::origin(); + let der = curve.der(0.5); + assert!(pt.cross(der)[2] > 0.0); + assert!(flags[1]); + flags[1] = false; + } + ShapesOpStatus::Or => { + let curve = bw[0].oriented_curve(); + let pt = curve.subs(0.5) - Point3::origin(); + let der = curve.der(0.5); + assert!(pt.cross(der)[2] < 0.0); + assert!(flags[2]); + flags[2] = false; + } + } + } + let mut flags = [true; 3]; + for i in 0..3 { + let bw = &geom_loops_store1[0][i]; + match bw.status() { + ShapesOpStatus::Unknown => { + let curve = bw[0].oriented_curve(); + assert_near!(curve.subs(0.0)[2], -1.0); + assert!(flags[0]); + flags[0] = false; + } + ShapesOpStatus::Or => { + let curve = bw[0].oriented_curve(); + let pt = curve.subs(0.5) - Point3::origin(); + let der = curve.der(0.5); + assert!(pt.cross(der)[2] < 0.0); + assert!(flags[1]); + flags[1] = false; + } + ShapesOpStatus::And => { + let curve = bw[0].oriented_curve(); + let pt = curve.subs(0.5) - Point3::origin(); + let der = curve.der(0.5); + assert!(pt.cross(der)[2] > 0.0); + assert!(flags[2]); + flags[2] = false; + } + } + } +} + +#[test] +fn rotated_intersection() { + // prepare geoetries + let arc00: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, 1.0, 0.0, 0.0), + Vector4::new(-1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc01: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(-1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, -1.0, 0.0, 0.0), + Vector4::new(1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc02: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, 0.0, -3.0, 1.0), + Vector4::new(-1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc10: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(0.0, -1.0, -1.0, 1.0), + Vector4::new(1.0, 0.0, 0.0, 0.0), + Vector4::new(0.0, 1.0, -1.0, 1.0), + ], + )) + .into(); + let arc11: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(0.0, 1.0, -1.0, 1.0), + Vector4::new(-1.0, 0.0, 0.0, 0.0), + Vector4::new(0.0, -1.0, -1.0, 1.0), + ], + )) + .into(); + let arc12: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(0.0, -1.0, -1.0, 1.0), + Vector4::new(0.0, 0.0, 3.0, 1.0), + Vector4::new(0.0, 1.0, -1.0, 1.0), + ], + )) + .into(); + let (surface0, surface1) = parabola_surfaces(); + + let v00 = Vertex::new(Point3::new(1.0, 0.0, 1.0)); + let v01 = Vertex::new(Point3::new(-1.0, 0.0, 1.0)); + let edge00 = Edge::new(&v00, &v01, arc00); + let edge01 = Edge::new(&v01, &v00, arc01); + let edge02 = Edge::new(&v00, &v01, arc02); + let wire00: Wire<_, _> = vec![edge00, edge02.inverse()].into(); + let wire01: Wire<_, _> = vec![edge01, edge02].into(); + let face00 = Face::new(vec![wire00], surface0.clone()); + let face01 = Face::new(vec![wire01], surface0); + let geom_shell0: Shell<_, _, _> = vec![face00.inverse(), face01.inverse()].into(); + + let v10 = Vertex::new(Point3::new(0.0, -1.0, -1.0)); + let v11 = Vertex::new(Point3::new(0.0, 1.0, -1.0)); + let edge10 = Edge::new(&v10, &v11, arc10); + let edge11 = Edge::new(&v11, &v10, arc11); + let edge12 = Edge::new(&v10, &v11, arc12); + let wire10: Wire<_, _> = vec![edge10, edge12.inverse()].into(); + let wire11: Wire<_, _> = vec![edge12, edge11].into(); + let face10 = Face::new(vec![wire10], surface1.clone()); + let face11 = Face::new(vec![wire11], surface1); + let geom_shell1: Shell<_, _, _> = vec![face10, face11].into(); + + let poly_shell0 = geom_shell0.triangulation(TOL); + let poly_shell1 = geom_shell1.triangulation(TOL); + let file = std::fs::File::create("polyshell0.obj").unwrap(); + obj::write(&poly_shell0.to_polygon(), file).unwrap(); + let file = std::fs::File::create("polyshell1.obj").unwrap(); + obj::write(&poly_shell1.to_polygon(), file).unwrap(); + let LoopsStoreQuadruple { + geom_loops_store0, + geom_loops_store1, + .. + } = create_loops_stores(&geom_shell0, &poly_shell0, &geom_shell1, &poly_shell1, TOL).unwrap(); + + let vertex_format = VertexDisplayFormat::AsPoint; + let edge_id_format = EdgeDisplayFormat::VerticesTupleAndID { vertex_format }; + let wire_id_format = WireDisplayFormat::EdgesListTuple { + edge_format: edge_id_format, + }; + let edge_geom_format = EdgeDisplayFormat::VerticesTupleAndCurve { vertex_format }; + let wire_geom_format = WireDisplayFormat::EdgesListTuple { + edge_format: edge_geom_format, + }; + assert_eq!(geom_loops_store0.len(), 2); + assert_eq!(geom_loops_store0[0].len(), 2); + assert!( + geom_loops_store0[0][0].is_closed(), + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert!( + geom_loops_store0[0][1].is_closed(), + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert!( + geom_loops_store0[0][0].is_geometric_consistent(), + "{:?}", + geom_loops_store0.display(wire_geom_format) + ); + assert!( + geom_loops_store0[0][1].is_geometric_consistent(), + "{:?}", + geom_loops_store0.display(wire_geom_format) + ); + let (a, b) = (geom_loops_store0[0][0].len(), geom_loops_store0[0][1].len()); + let compatible0 = match geom_loops_store0[0][0].status() { + ShapesOpStatus::And => a == 3, + ShapesOpStatus::Or => a == 5, + _ => false, + }; + let compatible1 = match geom_loops_store0[0][1].status() { + ShapesOpStatus::And => b == 3, + ShapesOpStatus::Or => b == 5, + _ => false, + }; + assert!( + compatible0 && compatible1 && a * b == 15, + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert_eq!(geom_loops_store0[1].len(), 2); + assert!( + geom_loops_store0[1][0].is_closed(), + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert!( + geom_loops_store0[1][1].is_closed(), + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert!( + geom_loops_store0[1][0].is_geometric_consistent(), + "{:?}", + geom_loops_store0.display(wire_geom_format) + ); + assert!( + geom_loops_store0[1][1].is_geometric_consistent(), + "{:?}", + geom_loops_store0.display(wire_geom_format) + ); + let (a, b) = (geom_loops_store0[1][0].len(), geom_loops_store0[1][1].len()); + let compatible0 = match geom_loops_store0[1][0].status() { + ShapesOpStatus::And => a == 3, + ShapesOpStatus::Or => a == 5, + _ => false, + }; + let compatible1 = match geom_loops_store0[1][1].status() { + ShapesOpStatus::And => b == 3, + ShapesOpStatus::Or => b == 5, + _ => false, + }; + assert!( + compatible0 && compatible1 && a * b == 15, + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert_eq!(geom_loops_store1.len(), 2); + assert_eq!(geom_loops_store1[0].len(), 2); + assert!( + geom_loops_store1[0][0].is_closed(), + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert!( + geom_loops_store1[0][1].is_closed(), + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert!( + geom_loops_store1[0][0].is_geometric_consistent(), + "{:?}", + geom_loops_store1.display(wire_geom_format) + ); + assert!( + geom_loops_store1[0][1].is_geometric_consistent(), + "{:?}", + geom_loops_store1.display(wire_geom_format) + ); + let (a, b) = (geom_loops_store1[0][0].len(), geom_loops_store1[0][1].len()); + let compatible0 = match geom_loops_store1[0][0].status() { + ShapesOpStatus::And => a == 3, + ShapesOpStatus::Or => a == 5, + _ => false, + }; + let compatible1 = match geom_loops_store1[0][1].status() { + ShapesOpStatus::And => b == 3, + ShapesOpStatus::Or => b == 5, + _ => false, + }; + assert!( + compatible0 && compatible1 && a * b == 15, + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert_eq!(geom_loops_store1[1].len(), 2); + assert!( + geom_loops_store1[1][0].is_closed(), + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert!( + geom_loops_store1[1][1].is_closed(), + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert!( + geom_loops_store1[1][0].is_geometric_consistent(), + "{:?}", + geom_loops_store1.display(wire_geom_format) + ); + assert!( + geom_loops_store1[1][1].is_geometric_consistent(), + "{:?}", + geom_loops_store1.display(wire_geom_format) + ); + let (a, b) = (geom_loops_store1[1][0].len(), geom_loops_store1[1][1].len()); + let compatible0 = match geom_loops_store1[1][0].status() { + ShapesOpStatus::And => a == 3, + ShapesOpStatus::Or => a == 5, + _ => false, + }; + let compatible1 = match geom_loops_store1[1][1].status() { + ShapesOpStatus::And => b == 3, + ShapesOpStatus::Or => b == 5, + _ => false, + }; + assert!( + compatible0 && compatible1 && a * b == 15, + "{:?}", + geom_loops_store1.display(wire_id_format) + ); +} + +#[test] +fn crossing_edges() { + let arc00: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, 1.0, 0.0, 0.0), + Vector4::new(-1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc01: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(-1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, -1.0, 0.0, 0.0), + Vector4::new(1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc02: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, 1.0, 1.0), + Vector4::new(0.0, 0.0, -3.0, 1.0), + Vector4::new(-1.0, 0.0, 1.0, 1.0), + ], + )) + .into(); + let arc10: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, -1.0, 1.0), + Vector4::new(0.0, 1.0, 0.0, 0.0), + Vector4::new(-1.0, 0.0, -1.0, 1.0), + ], + )) + .into(); + let arc11: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(-1.0, 0.0, -1.0, 1.0), + Vector4::new(0.0, -1.0, 0.0, 0.0), + Vector4::new(1.0, 0.0, -1.0, 1.0), + ], + )) + .into(); + let arc12: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new( + KnotVec::bezier_knot(2), + vec![ + Vector4::new(1.0, 0.0, -1.0, 1.0), + Vector4::new(0.0, 0.0, 3.0, 1.0), + Vector4::new(-1.0, 0.0, -1.0, 1.0), + ], + )) + .into(); + let (surface0, surface1) = parabola_surfaces(); + + let v00 = Vertex::new(Point3::new(1.0, 0.0, 1.0)); + let v01 = Vertex::new(Point3::new(-1.0, 0.0, 1.0)); + let edge00 = Edge::new(&v00, &v01, arc00); + let edge01 = Edge::new(&v01, &v00, arc01); + let edge02 = Edge::new(&v00, &v01, arc02); + let wire00: Wire<_, _> = vec![edge00, edge02.inverse()].into(); + let wire01: Wire<_, _> = vec![edge01, edge02].into(); + let face00 = Face::new(vec![wire00], surface0.clone()); + let face01 = Face::new(vec![wire01], surface0); + let geom_shell0: Shell<_, _, _> = vec![face00.inverse(), face01.inverse()].into(); + + let v10 = Vertex::new(Point3::new(1.0, 0.0, -1.0)); + let v11 = Vertex::new(Point3::new(-1.0, 0.0, -1.0)); + let edge10 = Edge::new(&v10, &v11, arc10); + let edge11 = Edge::new(&v11, &v10, arc11); + let edge12 = Edge::new(&v10, &v11, arc12); + let wire10: Wire<_, _> = vec![edge10, edge12.inverse()].into(); + let wire11: Wire<_, _> = vec![edge11, edge12].into(); + let face10 = Face::new(vec![wire10], surface1.clone()); + let face11 = Face::new(vec![wire11], surface1); + let geom_shell1: Shell<_, _, _> = vec![face10, face11].into(); + + let poly_shell0 = geom_shell0.triangulation(TOL); + let poly_shell1 = geom_shell1.triangulation(TOL); + let LoopsStoreQuadruple { + geom_loops_store0, + geom_loops_store1, + .. + } = create_loops_stores(&geom_shell0, &poly_shell0, &geom_shell1, &poly_shell1, TOL).unwrap(); + + let vertex_format = VertexDisplayFormat::AsPoint; + let edge_id_format = EdgeDisplayFormat::VerticesTupleAndID { vertex_format }; + let wire_id_format = WireDisplayFormat::EdgesListTuple { + edge_format: edge_id_format, + }; + let edge_geom_format = EdgeDisplayFormat::VerticesTupleAndCurve { vertex_format }; + let wire_geom_format = WireDisplayFormat::EdgesListTuple { + edge_format: edge_geom_format, + }; + assert_eq!(geom_loops_store0.len(), 2); + assert_eq!(geom_loops_store0[0].len(), 2); + assert!( + geom_loops_store0[0][0].is_closed(), + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert!( + geom_loops_store0[0][1].is_closed(), + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert!( + geom_loops_store0[0][0].is_geometric_consistent(), + "{:?}", + geom_loops_store0.display(wire_geom_format) + ); + assert!( + geom_loops_store0[0][1].is_geometric_consistent(), + "{:?}", + geom_loops_store0.display(wire_geom_format) + ); + let (a, b) = (geom_loops_store0[0][0].len(), geom_loops_store0[0][1].len()); + let compatible0 = match geom_loops_store0[0][0].status() { + ShapesOpStatus::And => a == 2, + ShapesOpStatus::Or => a == 4, + _ => false, + }; + let compatible1 = match geom_loops_store0[0][1].status() { + ShapesOpStatus::And => b == 2, + ShapesOpStatus::Or => b == 4, + _ => false, + }; + assert!( + compatible0 && compatible1 && a * b == 8, + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert_eq!(geom_loops_store0[1].len(), 2); + assert!( + geom_loops_store0[1][0].is_closed(), + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert!( + geom_loops_store0[1][1].is_closed(), + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert!( + geom_loops_store0[1][0].is_geometric_consistent(), + "{:?}", + geom_loops_store0.display(wire_geom_format) + ); + assert!( + geom_loops_store0[1][1].is_geometric_consistent(), + "{:?}", + geom_loops_store0.display(wire_geom_format) + ); + let (a, b) = (geom_loops_store0[1][0].len(), geom_loops_store0[1][1].len()); + let compatible0 = match geom_loops_store0[1][0].status() { + ShapesOpStatus::And => a == 2, + ShapesOpStatus::Or => a == 4, + _ => false, + }; + let compatible1 = match geom_loops_store0[1][1].status() { + ShapesOpStatus::And => b == 2, + ShapesOpStatus::Or => b == 4, + _ => false, + }; + assert!( + compatible0 && compatible1 && a * b == 8, + "{:?}", + geom_loops_store0.display(wire_id_format) + ); + assert_eq!(geom_loops_store1.len(), 2); + assert_eq!(geom_loops_store1[0].len(), 2); + assert!( + geom_loops_store1[0][0].is_closed(), + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert!( + geom_loops_store1[0][1].is_closed(), + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert!( + geom_loops_store1[0][0].is_geometric_consistent(), + "{:?}", + geom_loops_store1.display(wire_geom_format) + ); + assert!( + geom_loops_store1[0][1].is_geometric_consistent(), + "{:?}", + geom_loops_store1.display(wire_geom_format) + ); + let (a, b) = (geom_loops_store1[0][0].len(), geom_loops_store1[0][1].len()); + let compatible0 = match geom_loops_store1[0][0].status() { + ShapesOpStatus::And => a == 2, + ShapesOpStatus::Or => a == 4, + _ => false, + }; + let compatible1 = match geom_loops_store1[0][1].status() { + ShapesOpStatus::And => b == 2, + ShapesOpStatus::Or => b == 4, + _ => false, + }; + assert!( + compatible0 && compatible1 && a * b == 8, + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert_eq!(geom_loops_store1[1].len(), 2); + assert!( + geom_loops_store1[1][0].is_closed(), + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert!( + geom_loops_store1[1][1].is_closed(), + "{:?}", + geom_loops_store1.display(wire_id_format) + ); + assert!( + geom_loops_store1[1][0].is_geometric_consistent(), + "{:?}", + geom_loops_store1.display(wire_geom_format) + ); + assert!( + geom_loops_store1[1][1].is_geometric_consistent(), + "{:?}", + geom_loops_store1.display(wire_geom_format) + ); + let (a, b) = (geom_loops_store1[1][0].len(), geom_loops_store1[1][1].len()); + let compatible0 = match geom_loops_store1[1][0].status() { + ShapesOpStatus::And => a == 2, + ShapesOpStatus::Or => a == 4, + _ => false, + }; + let compatible1 = match geom_loops_store1[1][1].status() { + ShapesOpStatus::And => b == 2, + ShapesOpStatus::Or => b == 4, + _ => false, + }; + assert!( + compatible0 && compatible1 && a * b == 8, + "{:?}", + geom_loops_store0.display(wire_id_format) + ); +} diff --git a/vendor/truck-shapeops/src/transversal/mod.rs b/vendor/truck-shapeops/src/transversal/mod.rs new file mode 100644 index 00000000..de442179 --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/mod.rs @@ -0,0 +1,7 @@ +mod divide_face; +mod faces_classification; +mod integrate; +mod intersection_curve; +mod loops_store; +mod polyline_construction; +pub use integrate::{and, or, ShapeOpsCurve, ShapeOpsSurface}; diff --git a/vendor/truck-shapeops/src/transversal/polyline_construction/mod.rs b/vendor/truck-shapeops/src/transversal/polyline_construction/mod.rs new file mode 100644 index 00000000..c6862443 --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/polyline_construction/mod.rs @@ -0,0 +1,117 @@ +use rustc_hash::{FxHashMap as HashMap, FxHashSet as HashSet}; +use std::collections::VecDeque; +use truck_base::{cgmath64::*, tolerance::*}; +use truck_meshalgo::prelude::PolylineCurve; + +pub fn construct_polylines(lines: &[(Point3, Point3)]) -> Vec> { + let mut graph: Graph = lines.iter().collect(); + let mut res = Vec::new(); + while !graph.is_empty() { + let (mut idx, node) = graph.get_one(); + let mut wire: VecDeque<_> = vec![node.coord].into(); + while let Some((idx0, pt)) = graph.get_a_next_node(idx) { + idx = idx0; + wire.push_back(pt); + } + let mut idx = PointIndex::from(wire[0]); + while let Some((idx0, pt)) = graph.get_a_next_node(idx) { + idx = idx0; + wire.push_front(pt); + } + res.push(PolylineCurve(wire.into())); + } + res +} + +#[derive(Clone, Debug, Copy, Hash, PartialEq, Eq)] +struct PointIndex([i64; 3]); + +impl From for PointIndex { + #[inline(always)] + fn from(pt: Point3) -> PointIndex { + let idx = pt.add_element_wise(TOLERANCE) / (2.0 * TOLERANCE); + PointIndex(idx.cast::().unwrap().into()) + } +} + +struct Node { + coord: Point3, + adjacency: HashSet, +} + +impl Node { + #[inline(always)] + fn new(coord: Point3, adjacency: HashSet) -> Node { Node { coord, adjacency } } + + fn pop_one_adjacency(&mut self) -> PointIndex { + let idx = *self.adjacency.iter().next().unwrap(); + self.adjacency.remove(&idx); + idx + } +} + +struct Graph(HashMap); + +impl std::ops::Deref for Graph { + type Target = HashMap; + #[inline(always)] + fn deref(&self) -> &Self::Target { &self.0 } +} + +impl std::ops::DerefMut for Graph { + #[inline(always)] + fn deref_mut(&mut self) -> &mut Self::Target { &mut self.0 } +} + +impl Graph { + fn add_half_edge(&mut self, pt0: Point3, pt1: Point3) { + let idx0 = pt0.into(); + let idx1 = pt1.into(); + if let Some(node) = self.get_mut(&idx0) { + node.adjacency.insert(idx1); + } else { + let mut set = HashSet::default(); + set.insert(idx1); + self.insert(idx0, Node::new(pt0, set)); + } + } + + fn add_edge(&mut self, line: (Point3, Point3)) { + if !line.0.near(&line.1) { + self.add_half_edge(line.0, line.1); + self.add_half_edge(line.1, line.0); + } + } + + #[inline(always)] + fn get_one(&self) -> (PointIndex, &Node) { + let (idx, node) = self.iter().next().unwrap(); + (*idx, node) + } + + fn get_a_next_node(&mut self, idx: PointIndex) -> Option<(PointIndex, Point3)> { + let node = self.get_mut(&idx)?; + let idx0 = node.pop_one_adjacency(); + if node.adjacency.is_empty() { + self.remove(&idx); + } + let node = self.get_mut(&idx0)?; + node.adjacency.remove(&idx); + let pt = node.coord; + if node.adjacency.is_empty() { + self.remove(&idx0); + } + Some((idx0, pt)) + } +} + +impl<'a> FromIterator<&'a (Point3, Point3)> for Graph { + fn from_iter>(iter: I) -> Graph { + let mut res = Graph(HashMap::default()); + iter.into_iter().for_each(|line| res.add_edge(*line)); + res + } +} + +#[cfg(test)] +mod tests; diff --git a/vendor/truck-shapeops/src/transversal/polyline_construction/tests.rs b/vendor/truck-shapeops/src/transversal/polyline_construction/tests.rs new file mode 100644 index 00000000..4ec6fb08 --- /dev/null +++ b/vendor/truck-shapeops/src/transversal/polyline_construction/tests.rs @@ -0,0 +1,113 @@ +use super::*; + +#[test] +fn construct_polylines_positive0() { + let lines = vec![ + (Point3::new(0.0, 0.0, 0.0), Point3::new(1.0, 0.0, 0.0)), + (Point3::new(1.0, 0.0, 0.0), Point3::new(0.0, 1.0, 0.0)), + (Point3::new(1.0, 1.0, 0.0), Point3::new(0.0, 0.0, 1.0)), + (Point3::new(0.0, 1.0, 1.0), Point3::new(1.0, 1.0, 1.0)), + (Point3::new(0.0, 0.0, 1.0), Point3::new(1.0, 0.0, 1.0)), + (Point3::new(0.0, 1.0, 0.0), Point3::new(1.0, 1.0, 0.0)), + (Point3::new(1.0, 1.0, 1.0), Point3::new(0.0, 0.0, 0.0)), + (Point3::new(1.0, 0.0, 1.0), Point3::new(0.0, 1.0, 1.0)), + ]; + let polyline = construct_polylines(&lines); + assert_eq!(polyline.len(), 1); + assert_eq!(polyline[0].len(), 9); + + let mut sign = None; + for line in polyline[0].windows(2) { + let a = line[0][0] + line[0][1] * 2.0 + line[0][2] * 4.0; + let b = line[1][0] + line[1][1] * 2.0 + line[1][2] * 4.0; + let x = b - a; + assert!(f64::abs(x) == 1.0 || f64::abs(x) == 7.0); + let s = f64::signum(x * (x - 2.0) * (x + 2.0)); + if let Some(sign) = sign { + assert!(s == sign); + } else { + sign = Some(s); + } + } +} + +#[test] +fn construct_polylines_positive1() { + let lines = vec![ + (Point3::new(0.0, 0.0, 0.0), Point3::new(1.0, 0.0, 0.0)), + (Point3::new(1.0, 0.0, 0.0), Point3::new(1.0, 1.0, 0.0)), + (Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 1.0, 1.0)), + (Point3::new(1.0, 1.0, 1.0), Point3::new(0.0, 1.0, 1.0)), + (Point3::new(1.0, 1.0, 0.0), Point3::new(0.0, 1.0, 0.0)), + (Point3::new(0.0, 1.0, 0.0), Point3::new(0.0, 0.0, 0.0)), + (Point3::new(0.0, 0.0, 1.0), Point3::new(1.0, 0.0, 1.0)), + (Point3::new(0.0, 1.0, 1.0), Point3::new(0.0, 0.0, 1.0)), + ]; + let polyline = construct_polylines(&lines); + assert_eq!(polyline.len(), 2); + assert_eq!(polyline[0].len(), 5); + assert_eq!(polyline[1].len(), 5); +} + +#[test] +fn construct_polylines_positive2() { + let lines = vec![ + (Point3::new(0.0, 0.0, 0.0), Point3::new(1.0, 0.0, 0.0)), + (Point3::new(1.0, 0.0, 0.0), Point3::new(0.0, 1.0, 0.0)), + (Point3::new(1.0, 1.0, 0.0), Point3::new(0.0, 0.0, 1.0)), + (Point3::new(0.0, 1.0, 1.0), Point3::new(1.0, 1.0, 1.0)), + (Point3::new(0.0, 0.0, 1.0), Point3::new(1.0, 0.0, 1.0)), + (Point3::new(1.0, 1.0, 0.0), Point3::new(1.0, 1.0, 0.0)), + (Point3::new(0.0, 1.0, 0.0), Point3::new(1.0, 1.0, 0.0)), + (Point3::new(1.0, 1.0, 1.0), Point3::new(0.0, 0.0, 0.0)), + (Point3::new(1.0, 0.0, 1.0), Point3::new(0.0, 1.0, 1.0)), + ]; + let polyline = construct_polylines(&lines); + assert_eq!(polyline.len(), 1); + assert_eq!(polyline[0].len(), 9); + + let mut sign = None; + for line in polyline[0].windows(2) { + let a = line[0][0] + line[0][1] * 2.0 + line[0][2] * 4.0; + let b = line[1][0] + line[1][1] * 2.0 + line[1][2] * 4.0; + let x = b - a; + assert!(f64::abs(x) == 1.0 || f64::abs(x) == 7.0); + let s = f64::signum(x * (x - 2.0) * (x + 2.0)); + if let Some(sign) = sign { + assert!(s == sign); + } else { + sign = Some(s); + } + } +} + +#[test] +fn construct_polylines_positive3() { + let lines = vec![ + (Point3::new(0.0, 0.0, 0.0), Point3::new(1.0, 0.0, 0.0)), + (Point3::new(1.0, 0.0, 0.0), Point3::new(0.0, 1.0, 0.0)), + (Point3::new(1.0, 1.0, 0.0), Point3::new(0.0, 0.0, 1.0)), + (Point3::new(0.0, 1.0, 1.0), Point3::new(1.0, 1.0, 1.0)), + (Point3::new(0.0, 0.0, 1.0), Point3::new(1.0, 0.0, 1.0)), + (Point3::new(1.0, 1.0, 0.0), Point3::new(1.0, 1.0, 0.0)), + (Point3::new(0.0, 1.0, 0.0), Point3::new(1.0, 1.0, 0.0)), + (Point3::new(1.0, 0.0, 1.0), Point3::new(0.0, 1.0, 1.0)), + ]; + let polyline = construct_polylines(&lines); + assert_eq!(polyline.len(), 1); + assert_eq!(polyline[0].len(), 8); + + let mut sign = None; + for line in polyline[0].windows(2) { + let a = line[0][0] + line[0][1] * 2.0 + line[0][2] * 4.0; + let b = line[1][0] + line[1][1] * 2.0 + line[1][2] * 4.0; + let x = b - a; + assert!(f64::abs(x) == 1.0); + let s = f64::signum(x * (x - 2.0) * (x + 2.0)); + if let Some(sign) = sign { + assert!(s == sign); + } else { + sign = Some(s); + } + } +}