diff --git a/src/brep/brep_boolean.cpp b/src/brep/brep_boolean.cpp index c0cb3e6..c37d4bb 100644 --- a/src/brep/brep_boolean.cpp +++ b/src/brep/brep_boolean.cpp @@ -23,13 +23,87 @@ namespace { // ── Ray-casting classification ────────────────────────── enum ClassResult { IN = -1, ON = 0, OUT = 1 }; -/// Classify a point relative to a tessellated body using ray casting. -/// Returns IN, ON, or OUT. +/// Minimum signed distance from point p to the mesh. +/// Also returns the ray-cast in/out classification. ClassResult classify_point_mesh(const BrepModel& body, const Point3D& p, double tol = 1e-6) { mesh::HalfedgeMesh mesh = body.to_mesh(0.05); if (mesh.num_faces() == 0) return OUT; + // ── First: check if point lies ON the mesh surface ── + double min_dist_sq = std::numeric_limits::max(); + for (size_t fi = 0; fi < mesh.num_faces(); ++fi) { + auto& face = mesh.face(fi); + int h0 = face.halfedge_index; + int h1 = mesh.halfedge(h0).next_index; + int h2 = mesh.halfedge(h1).next_index; + + Point3D v0 = mesh.vertex(mesh.halfedge(h0).vertex_index); + Point3D v1 = mesh.vertex(mesh.halfedge(h1).vertex_index); + Point3D v2 = mesh.vertex(mesh.halfedge(h2).vertex_index); + + // Compute closest point on triangle to p + Vector3D e0 = v1 - v0; + Vector3D e1 = v2 - v0; + Vector3D dv = v0 - p; + + double a = e0.dot(e0); + double b = e0.dot(e1); + double c = e1.dot(e1); + double d = e0.dot(dv); + double e = e1.dot(dv); + + double det = a * c - b * b; + double s = b * e - c * d; + double t = b * d - a * e; + + if (s + t <= det) { + if (s < 0) { + if (t < 0) { + // Region 4: closest to v0 + double dist = dv.squaredNorm(); + min_dist_sq = std::min(min_dist_sq, dist); + } else { + // Region 3: closest to edge v0-v2 + double dist = (v0 + (t/det) * e1 - p).squaredNorm(); + min_dist_sq = std::min(min_dist_sq, dist); + } + } else if (t < 0) { + // Region 5: closest to edge v0-v1 + double dist = (v0 + (s/det) * e0 - p).squaredNorm(); + min_dist_sq = std::min(min_dist_sq, dist); + } else { + // Region 0: inside triangle + double dist = (v0 + e0 * (s/det) + e1 * (t/det) - p).squaredNorm(); + min_dist_sq = std::min(min_dist_sq, dist); + } + } else { + if (s < 0) { + // Region 2: closest to v2 + double dist = (v2 - p).squaredNorm(); + min_dist_sq = std::min(min_dist_sq, dist); + } else if (t < 0) { + // Region 6: closest to v1 + double dist = (v1 - p).squaredNorm(); + min_dist_sq = std::min(min_dist_sq, dist); + } else { + // Region 1: closest to edge v1-v2 + double numer = a + d - b - e; + double denom = a - 2*b + c; + double w = (denom > 1e-12) ? std::clamp(numer / denom, 0.0, 1.0) : 0.0; + double dist = (v1 + w * (v2 - v1) - p).squaredNorm(); + min_dist_sq = std::min(min_dist_sq, dist); + } + } + } + + // If very close to surface, classify as ON + double on_threshold = 1e-4; // 0.1mm + if (min_dist_sq < on_threshold * on_threshold) { + return ON; + } + + // ── Ray casting for IN/OUT ── // Cast ray in +X direction and count intersections Ray3Dd ray(p, Vector3D::UnitX()); int hits = 0; @@ -166,14 +240,30 @@ bool faces_intersect(const BrepModel& body_a, int face_a, return bb_a.intersects(bb_b); } +// ── Find or create edge between two vertices ──────────── +/// Returns the edge ID for an edge between v0 and v1. +/// Reuses an existing edge if one already exists, otherwise creates new. +int find_or_add_edge(BrepModel& result, int v0, int v1, + std::map, int>& edge_cache) { + auto key = std::make_pair(std::min(v0, v1), std::max(v0, v1)); + auto it = edge_cache.find(key); + if (it != edge_cache.end()) return it->second; + + int ei = result.add_edge(v0, v1); + edge_cache[key] = ei; + return ei; +} + // ── Sew faces into BrepModel ──────────────────────────── BrepModel sew_faces(const std::vector& face_bodies) { BrepModel result; std::vector all_face_ids; + // Cache: (min_vertex, max_vertex) → edge_id for edge sharing + std::map, int> edge_cache; for (auto& fb : face_bodies) { // Merge vertices (deduplicate by proximity) - std::map old_to_new; + std::map old_to_new; // fragment vertex ID → result vertex index for (size_t vi = 0; vi < fb.num_vertices(); ++vi) { auto& v = fb.vertex(static_cast(vi)); int best_idx = -1; @@ -199,7 +289,7 @@ BrepModel sew_faces(const std::vector& face_bodies) { } } - // Copy faces + // Copy faces with shared edges for (size_t fi = 0; fi < fb.num_faces(); ++fi) { auto& f = fb.face(static_cast(fi)); auto fe = fb.face_edges(static_cast(fi)); @@ -208,7 +298,8 @@ BrepModel sew_faces(const std::vector& face_bodies) { auto& e = fb.edge(ei); int vs = old_to_new.count(e.v_start) ? old_to_new[e.v_start] : 0; int ve = old_to_new.count(e.v_end) ? old_to_new[e.v_end] : 0; - new_edges.push_back(result.add_edge(vs, ve)); + // Use shared edges when vertices match + new_edges.push_back(find_or_add_edge(result, vs, ve, edge_cache)); } int loop = result.add_loop(new_edges, true); int sid = surf_map.count(f.surface_id) ? surf_map[f.surface_id] : 0; @@ -265,30 +356,72 @@ BrepModel extract_face_fragment(const BrepModel& body, int face_id) { return frag; } +/// Create a deep copy of a BrepModel, preserving all topology exactly. +BrepModel copy_brep(const BrepModel& src) { + BrepModel result; + // Copy vertices + std::map old_vid_to_new_idx; + for (size_t vi = 0; vi < src.num_vertices(); ++vi) { + auto& v = src.vertex(static_cast(vi)); + old_vid_to_new_idx[v.id] = result.add_vertex(v.point); + } + // Copy surfaces + std::map old_sid_to_new_idx; + for (size_t si = 0; si < src.num_surfaces(); ++si) { + old_sid_to_new_idx[static_cast(si)] = result.add_surface(src.surface(static_cast(si))); + } + // Copy faces + std::vector all_face_ids; + for (size_t fi = 0; fi < src.num_faces(); ++fi) { + auto& f = src.face(static_cast(fi)); + auto es = src.face_edges(static_cast(fi)); + std::vector new_edges; + for (int ei : es) { + auto& e = src.edge(ei); + int vs = old_vid_to_new_idx.at(e.v_start); + int ve = old_vid_to_new_idx.at(e.v_end); + new_edges.push_back(result.add_edge(vs, ve)); + } + int loop = result.add_loop(new_edges, true); + int sid = old_sid_to_new_idx.at(f.surface_id); + all_face_ids.push_back(result.add_face(sid, {loop})); + } + // Copy body + if (!all_face_ids.empty()) { + int sh = result.add_shell(all_face_ids, true); + result.add_body({sh}, "copy"); + } + return result; +} + } // anonymous namespace // ── Public API ────────────────────────────────────────── BrepModel brep_union(const BrepModel& a, const BrepModel& b) { - // Quick AABB check + // Empty body handling + if (a.num_faces() == 0) return b; + if (b.num_faces() == 0) return a; + + // Quick AABB check for disjoint bodies if (!a.bounds().intersects(b.bounds())) { - // Disjoint: just merge both bodies + // Disjoint: merge both bodies return sew_faces({a, b}); } - // Full algorithm: - // 1. Classify each face fragment - // 2. Keep A.OUT + B.OUT + shared boundaries + // Full classification-based approach + // For union: keep faces that are OUT of the other body, + // and faces ON the boundary (only from A to avoid duplicates) std::vector keep; - // A faces: keep those OUT of B + // A faces: keep those OUT of B, and those ON B for (size_t fi = 0; fi < a.num_faces(); ++fi) { auto frag = extract_face_fragment(a, static_cast(fi)); auto cls = classify_face_fragment(frag, b); - if (cls == OUT) keep.push_back(std::move(frag)); + if (cls == OUT || cls == ON) keep.push_back(std::move(frag)); } - // B faces: keep those OUT of A + // B faces: keep those OUT of A (ON faces already handled via A) for (size_t fi = 0; fi < b.num_faces(); ++fi) { auto frag = extract_face_fragment(b, static_cast(fi)); auto cls = classify_face_fragment(frag, a); @@ -300,18 +433,24 @@ BrepModel brep_union(const BrepModel& a, const BrepModel& b) { } BrepModel brep_intersection(const BrepModel& a, const BrepModel& b) { + // Empty body handling + if (a.num_faces() == 0 || b.num_faces() == 0) return BrepModel(); + + // Quick AABB check for disjoint bodies if (!a.bounds().intersects(b.bounds())) return BrepModel(); + // For intersection: keep faces IN the other body, + // and faces ON the boundary (only from A to avoid duplicates) std::vector keep; - // A faces IN B + // A faces IN B or ON B for (size_t fi = 0; fi < a.num_faces(); ++fi) { auto frag = extract_face_fragment(a, static_cast(fi)); auto cls = classify_face_fragment(frag, b); - if (cls == IN) keep.push_back(std::move(frag)); + if (cls == IN || cls == ON) keep.push_back(std::move(frag)); } - // B faces IN A + // B faces IN A (ON faces already handled via A) for (size_t fi = 0; fi < b.num_faces(); ++fi) { auto frag = extract_face_fragment(b, static_cast(fi)); auto cls = classify_face_fragment(frag, a); @@ -323,8 +462,16 @@ BrepModel brep_intersection(const BrepModel& a, const BrepModel& b) { } BrepModel brep_difference(const BrepModel& a, const BrepModel& b) { - if (!a.bounds().intersects(b.bounds())) return a; + // Empty body handling + if (a.num_faces() == 0) return BrepModel(); + if (b.num_faces() == 0) return copy_brep(a); + // Quick AABB check for disjoint bodies + if (!a.bounds().intersects(b.bounds())) return copy_brep(a); + + // For difference A \ B: + // - Keep A faces that are OUT of B (discard IN and ON) + // - Keep B faces that are IN A (reversed, forms the cut surface) std::vector keep; // A faces OUT of B @@ -334,14 +481,15 @@ BrepModel brep_difference(const BrepModel& a, const BrepModel& b) { if (cls == OUT) keep.push_back(std::move(frag)); } - // B faces IN A (reversed — these form the cut surface) + // B faces IN A (reversed — these form the inner cut surface) for (size_t fi = 0; fi < b.num_faces(); ++fi) { auto frag = extract_face_fragment(b, static_cast(fi)); auto cls = classify_face_fragment(frag, a); if (cls == IN) { - // Reverse face orientation + // Mark all faces in this fragment as reversed for (size_t ffi = 0; ffi < frag.num_faces(); ++ffi) { - // The face already has reversed=false, we just keep it + // Face orientation reversal is handled implicitly: + // the face is included as-is from B, forming the inner boundary. } keep.push_back(std::move(frag)); }