diff --git a/src/cpp/web-ifc/geometry/IfcGeometryProcessor.cpp b/src/cpp/web-ifc/geometry/IfcGeometryProcessor.cpp index c28eb18d..f17c2d2a 100644 --- a/src/cpp/web-ifc/geometry/IfcGeometryProcessor.cpp +++ b/src/cpp/web-ifc/geometry/IfcGeometryProcessor.cpp @@ -215,6 +215,10 @@ namespace webifc::geometry std::vector geomVector = {geom}; // Wrap 'geom' in a vector fusedVoids = BoolProcess(std::vector{fusedVoids}, geomVector, "UNION", _settings); } + if (fusedVoids.numFaces > _settings._CSG_MAX_NUM_FACES) { + spdlog::warn("High number of faces ({}) in voids for element {}, skipping further CSG operations", fusedVoids.numFaces, expressID); + break; + } } joinedVoidGeoms.push_back(fusedVoids); @@ -339,6 +343,39 @@ namespace webifc::geometry return mesh; } + size_t numFacesGeoms = 0; + for (auto& geom : flatFirstMeshes) { + numFacesGeoms += geom.numFaces; + } + for (auto& geom : flatSecondMeshes) { + numFacesGeoms += geom.numFaces; + } + + if (numFacesGeoms > _settings._CSG_MAX_NUM_FACES) { + spdlog::warn("High number of faces in CSG operand ({}), skipping further CSG operations", numFacesGeoms); + + // Flatten the first operand (transformed and normalized, like in the normal path) + auto flatFirstMeshes = flatten(firstMesh, _expressIDToGeometry, normalizeMat); + + // Merge into a single geometry (equivalent to the first operand without boolean) + IfcGeometry resultMesh; + for (const auto& g : flatFirstMeshes) { + resultMesh.MergeGeometry(g); + } + + _expressIDToGeometry[expressID] = resultMesh; + mesh.hasGeometry = true; + mesh.transformation = glm::translate(origin); + + if (!mesh.hasColor && firstMesh.hasColor) { + mesh.hasColor = true; + mesh.color = firstMesh.color; + } + + // No need to set mesh = firstMesh; we're building a flat geometry here for consistency + return mesh; + } + IfcGeometry resultMesh = BoolProcess(flatFirstMeshes, flatSecondMeshes, std::string(op), _settings); _expressIDToGeometry[expressID] = resultMesh; diff --git a/src/cpp/web-ifc/geometry/IfcGeometryProcessor.h b/src/cpp/web-ifc/geometry/IfcGeometryProcessor.h index 1ef76932..7d5a7a45 100644 --- a/src/cpp/web-ifc/geometry/IfcGeometryProcessor.h +++ b/src/cpp/web-ifc/geometry/IfcGeometryProcessor.h @@ -34,6 +34,7 @@ namespace webifc::geometry double TOLERANCE_INSIDE_OUTSIDE_PERIMETER = 1.0E-10; double TOLERANCE_BOUNDING_BOX = 1.0E-02; uint16_t _BOOLEAN_UNION_THRESHOLD = 150; + uint16_t _CSG_MAX_NUM_FACES = 20000; }; class booleanManager diff --git a/src/cpp/web-ifc/geometry/operations/boolean-utils/shared-position.h b/src/cpp/web-ifc/geometry/operations/boolean-utils/shared-position.h index 32719ecc..2f2f2446 100644 --- a/src/cpp/web-ifc/geometry/operations/boolean-utils/shared-position.h +++ b/src/cpp/web-ifc/geometry/operations/boolean-utils/shared-position.h @@ -1,18 +1,14 @@ #pragma once - #include #include #include #include #include - #include - #pragma warning(push) #pragma warning(disable : 4267) #include #pragma warning(pop) - #include "geometry.h" #include "aabb.h" #include "bvh.h" @@ -23,35 +19,37 @@ #include "obj-exporter.h" #include "is-inside-mesh.h" #include "is-inside-boundary.h" - #define GLM_ENABLE_EXPERIMENTAL #include - using Vec2 = glm::dvec2; using Vec3 = glm::dvec3; - +namespace std { + template <> + struct hash { + size_t operator()(const glm::ivec3& v) const { + size_t h1 = std::hash()(v.x); + size_t h2 = std::hash()(v.y); + size_t h3 = std::hash()(v.z); + return h1 ^ (h2 << 1) ^ (h3 << 2); // Simple combine + } + }; +} namespace fuzzybools { - //============================================================================================ - struct PlaneBasis { Vec3 origin; - Vec3 up; Vec3 left; Vec3 right; - - Vec2 project(const Vec3 &pt) + Vec2 project(const Vec3& pt) { auto relative = pt - origin; return Vec2(glm::dot(relative, left), glm::dot(relative, right)); } }; - //============================================================================================ - struct ReferencePlane { size_t planeID; @@ -59,21 +57,18 @@ namespace fuzzybools size_t lineID; Vec2 location; }; - //============================================================================================ - struct ReferenceLine { size_t lineID; size_t pointID; double location; }; - //============================================================================================ //============================================================================================ /* - The vector direction is not necessarily a unit vector! It is the vector pointing - from the start of the line (identified with origin) to the end. Thus, + The vector direction is not necessarily a unit vector! It is the vector pointing + from the start of the line (identified with origin) to the end. Thus, end = origin + direction */ struct Line @@ -82,23 +77,18 @@ namespace fuzzybools size_t globalID; Vec3 origin; Vec3 direction; - //============================================================================================ - Line() { static size_t idcounter = 0; idcounter++; globalID = idcounter; } - //============================================================================================ - /* Does point a lie on the line defined by vectors origin and direction? - */ - bool IsPointOnLine(const Vec3 &a) const + bool IsPointOnLine(const Vec3& a) const { /* Vector d is the unit vector pointing along the line from origin towards @@ -107,7 +97,7 @@ namespace fuzzybools Vec3 d = glm::normalize(direction); Vec3 v = a - origin; /* - Drop a perpendicular from point a to the line. The quantity t + Drop a perpendicular from point a to the line. The quantity t is the distance from origin along the line to the foot of the perpendicular, denoted by p. */ @@ -119,43 +109,35 @@ namespace fuzzybools */ return glm::distance(p, a) < tolerancePointOnLine; } - //============================================================================================ - /* The original version of GetPosOnline seemed to assume that direction is a unit vector. This version normalises direction first. */ - - double GetPosOnLine(const Vec3 &pos) const + double GetPosOnLine(const Vec3& pos) const { Vec3 unitDirection = glm::normalize(direction); return glm::dot(pos - origin, unitDirection); } - Vec3 GetPosOnLine(const double dist) const { Vec3 unitDirection = glm::normalize(direction); return origin + dist * unitDirection; } - //============================================================================================ - - bool IsCollinear(const Line &other) const + bool IsCollinear(const Line& other) const { Vec3 unitDirection = glm::normalize(direction); Vec3 unitOtherDirection = glm::normalize(other.direction); return (equals(unitOtherDirection, unitDirection, toleranceCollinear) || equals(unitOtherDirection, -unitDirection, toleranceCollinear)); - // return (equals(other.direction, direction, toleranceCollinear) || equals(other.direction, -direction, toleranceCollinear)); + // return (equals(other.direction, direction, toleranceCollinear) || equals(other.direction, -direction, toleranceCollinear)); } - //============================================================================================ - /* - The original version of IsEqualTo compared dir and direction. This version compares + The original version of IsEqualTo compared dir and direction. This version compares normalised copies of dir and direction. */ - bool IsEqualTo(const Vec3 &pos, const Vec3 &dir) const + bool IsEqualTo(const Vec3& pos, const Vec3& dir) const { // check dir Vec3 unitDir = glm::normalize(dir); @@ -164,87 +146,79 @@ namespace fuzzybools { return false; } - // check pos if (!IsPointOnLine(pos)) { return false; } - return true; } - //============================================================================================ - void AddPointToLine(double dist, size_t id) { // check existing - for (auto &p : points) + for (auto& p : points) { if (p.second == id) return; } - // add new point points.push_back(std::make_pair(dist, id)); - // re-sort all std::sort( points.begin(), points.end(), - [&](const std::pair &left, const std::pair &right) + [&](const std::pair& left, const std::pair& right) { return left.first < right.first; }); } - //============================================================================================ - auto GetSegments() const { const auto makeSegments = [&](int i) - { - return std::make_pair(points[i - 1].second, points[i].second); - }; - + { + return std::make_pair(points[i - 1].second, points[i].second); + }; return std::views::iota(size_t(1), points.size()) | std::views::transform(makeSegments); } - + //============================================================================================ + bool ContainsPos(double pos, double tol = tolerancePointOnLine) const + { + if (points.empty()) return false; + for (size_t i = 1; i < points.size(); ++i) { + double minv = points[i - 1].first; + double maxv = points[i].first; + if (pos >= minv - tol && pos <= maxv + tol) return true; + } + return false; + } std::vector> points; - std::vector planes; }; - //============================================================================================ //============================================================================================ - struct Point { size_t id; size_t globalID; Vec3 location3D; - Point() { static size_t idcounter = 0; idcounter++; globalID = idcounter; } - //============================================================================================ - - bool operator==(const Vec3 &pt) + bool operator==(const Vec3& pt) { return equals(location3D, pt, toleranceVectorEquality); } - std::vector lines; std::vector planes; }; - //============================================================================================ //============================================================================================ - struct Plane { int refPlane = -1; @@ -252,64 +226,48 @@ namespace fuzzybools size_t globalID; double distance; Vec3 normal; - std::vector lines; AABB aabb; - //============================================================================================ - - void AddPoint(const Vec3 &pt) + void AddPoint(const Vec3& pt) { aabb.merge(pt); } - //============================================================================================ - Plane() { static size_t idcounter = 0; idcounter++; globalID = idcounter; } - //============================================================================================ - double round(double input) { input = std::fabs(input) < EPS_BIG ? 0.0 : input; input = std::fabs(input) < (1.0 - EPS_BIG) ? input : input > 0.0 ? 1.0 - : -1.0; + : -1.0; return input; } - //============================================================================================ - Vec3 round(Vec3 in) { in.x = round(in.x); in.y = round(in.y); in.z = round(in.z); - return in; } - //============================================================================================ - Vec3 GetDirection(Vec3 a, Vec3 b) { auto dir = b - a; return glm::normalize(dir); } - //============================================================================================ - - std::pair AddLine(const Point &a, const Point &b) + std::pair AddLine(const Point& a, const Point& b) { Vec3 pos = a.location3D; Vec3 dir = GetDirection(pos, b.location3D); - auto lineId = AddLine(pos, dir); - if (!lines[lineId.first].IsPointOnLine(a.location3D)) { if (messages) @@ -324,7 +282,6 @@ namespace fuzzybools printf("bad point in AddLine\n"); } } - if (!aabb.contains(a.location3D)) { if (messages) @@ -339,70 +296,52 @@ namespace fuzzybools printf("bad points in AddLine\n"); } } - lines[lineId.first].AddPointToLine(lines[lineId.first].GetPosOnLine(a.location3D), a.id); lines[lineId.first].AddPointToLine(lines[lineId.first].GetPosOnLine(b.location3D), b.id); - return lineId; } - //============================================================================================ - - std::pair AddLine(const Vec3 &pos, const Vec3 &dir) + std::pair AddLine(const Vec3& pos, const Vec3& dir) { - for (auto &line : lines) + for (auto& line : lines) { if (line.IsEqualTo(pos, dir)) { - return {line.id, false}; + return { line.id, false }; } } - Line l; l.id = lines.size(); l.origin = pos; - // l.direction = dir; Vec3 temp = glm::normalize(dir); l.direction = temp; - lines.push_back(l); - - return {l.id, true}; + return { l.id, true }; } - //============================================================================================ - void RemoveLastLine() { lines.pop_back(); } - //============================================================================================ - - bool IsEqualTo(const Vec3 &n, double d) + bool IsEqualTo(const Vec3& n, double d) { return (equals(normal, n, toleranceVectorEquality) && equals(distance, d, TOLERANCE_SCALAR_EQUALITY)); } - //============================================================================================ - - glm::dvec2 GetPosOnPlane(const glm::dvec3 &pos) + glm::dvec2 GetPosOnPlane(const glm::dvec3& pos) { return {}; } - //============================================================================================ - - bool HasOverlap(const std::pair &A, const std::pair &B) + bool HasOverlap(const std::pair& A, const std::pair& B) { return (A.first == B.first || A.first == B.second || A.second == B.first || A.second == B.second); } - //============================================================================================ - - void PutPointOnLines(Point &p) + void PutPointOnLines(Point& p) { - for (auto &l : lines) + for (auto& l : lines) { if (l.IsPointOnLine(p.location3D)) { @@ -414,89 +353,67 @@ namespace fuzzybools } } } - //============================================================================================ - /* Normal is assumed to be normalised. */ - bool IsPointOnPlane(const glm::dvec3 &pos) + bool IsPointOnPlane(const glm::dvec3& pos) { double d = glm::dot(normal, pos); double posLength = pos.length(); - // return equals(distance, d, toleranceVectorEquality * posLength); + // return equals(distance, d, toleranceVectorEquality * posLength); return equals(distance, d, toleranceVectorEquality); } - //============================================================================================ - PlaneBasis MakeBasis() { glm::dvec3 origin = normal * distance; glm::dvec3 up = normal; - glm::dvec3 worldUp = glm::dvec3(0, 1, 0); glm::dvec3 worldRight = glm::dvec3(1, 0, 0); - bool normalIsUp = equals(up, worldUp, EPS_SMALL) || equals(-up, worldUp, EPS_SMALL); glm::dvec3 left = normalIsUp ? glm::cross(up, worldRight) : glm::cross(up, worldUp); glm::dvec3 right = glm::cross(left, up); - PlaneBasis basis; - basis.origin = origin; basis.up = glm::normalize(up); basis.left = glm::normalize(left); basis.right = glm::normalize(right); - return basis; } }; - //============================================================================================ //============================================================================================ - struct Triangle { size_t id; - size_t a; size_t b; size_t c; - //============================================================================================ - void Flip() { auto temp = a; a = b; b = temp; } - //============================================================================================ - bool HasPoint(size_t p) { return a == p || b == p || c == p; } - //============================================================================================ - - bool IsNeighbour(Triangle &t) + bool IsNeighbour(Triangle& t) { return HasPoint(t.a) || HasPoint(t.b) || HasPoint(t.c); } - //============================================================================================ - - bool SamePoints(Triangle &other) + bool SamePoints(Triangle& other) { return other.HasPoint(a) && other.HasPoint(b) && other.HasPoint(c); } - //============================================================================================ - - size_t GetNotShared(Triangle &other) + size_t GetNotShared(Triangle& other) { if (!other.HasPoint(a)) { @@ -510,14 +427,11 @@ namespace fuzzybools { return c; } - // throw std::exception("Missing common point in GetNotShared"); } }; - //============================================================================================ //============================================================================================ - struct SegmentSet { std::vector> segments; @@ -525,12 +439,9 @@ namespace fuzzybools std::map, size_t> segmentCounts; std::vector irrelevantFaces; std::vector irrelevantFaces_toTest; - std::map>> planeSegments; std::map, size_t>> planeSegmentCounts; - //============================================================================================ - void AddSegment(size_t planeId, size_t a, size_t b) { if (a == b) @@ -541,73 +452,55 @@ namespace fuzzybools } return; } - auto seg = a < b ? std::make_pair(a, b) : std::make_pair(b, a); - segments.emplace_back(seg); segmentCounts[seg]++; - planeSegments[planeId].emplace_back(seg); planeSegmentCounts[planeId][seg]++; } - //============================================================================================ - void AddFace(size_t planeId, size_t a, size_t b, size_t c) { AddSegment(planeId, a, b); AddSegment(planeId, b, c); AddSegment(planeId, c, a); - Triangle t; t.id = triangles.size(); t.a = a; t.b = b; t.c = c; - triangles.push_back(t); } - //============================================================================================ - std::vector GetTrianglesWithPoint(size_t p) { std::vector returnTriangles; - - for (auto &t : triangles) + for (auto& t : triangles) { if (t.HasPoint(p)) { returnTriangles.push_back(t.id); } } - return returnTriangles; } - //============================================================================================ - std::vector GetTrianglesWithEdge(size_t a, size_t b) { std::vector returnTriangles; - - for (auto &t : triangles) + for (auto& t : triangles) { if (t.HasPoint(a) && t.HasPoint(b)) { returnTriangles.push_back(t.id); } } - return returnTriangles; } - //============================================================================================ - - std::vector, std::vector>> GetNeighbourTriangles(Triangle &triangle) + std::vector, std::vector>> GetNeighbourTriangles(Triangle& triangle) { std::vector, std::vector>> returnTriangles; - { auto tris = GetTrianglesWithEdge(triangle.a, triangle.b); returnTriangles.emplace_back(std::make_pair(triangle.a, triangle.b), tris); @@ -620,32 +513,26 @@ namespace fuzzybools auto tris = GetTrianglesWithEdge(triangle.c, triangle.a); returnTriangles.emplace_back(std::make_pair(triangle.c, triangle.a), tris); } - return returnTriangles; } - bool IsManifold() { std::vector> contours; - - for (auto &[pair, count] : segmentCounts) + for (auto& [pair, count] : segmentCounts) { if (count != 2) { contours.push_back(pair); } } - return contours.empty(); } - std::map>> GetContourSegments() { std::map>> contours; - - for (auto &[plane, segmentCounts] : planeSegmentCounts) + for (auto& [plane, segmentCounts] : planeSegmentCounts) { - for (auto &[pair, count] : segmentCounts) + for (auto& [pair, count] : segmentCounts) { if (count == 1) { @@ -653,61 +540,44 @@ namespace fuzzybools } } } - return contours; } }; /////////////////////////////////////////////////// - //============================================================================================ //============================================================================================ - struct SharedPosition { std::vector points; std::vector planes; - SegmentSet A; SegmentSet B; - // TODO: design flaw - const Geometry *_linkedA; - const Geometry *_linkedB; - + const Geometry* _linkedA; + const Geometry* _linkedB; Geometry relevantA; Geometry relevantB; - BVH relevantBVHA; BVH relevantBVHB; - //============================================================================================ - - // assumes all triangleIds are connected to base with an edge and are flipped correctly - size_t FindUppermostTriangleId(Triangle &base, const std::vector &triangleIds) + // assumes all triangleIds are connected to base with an edge and are flipped correctly + size_t FindUppermostTriangleId(Triangle& base, const std::vector& triangleIds) { if (triangleIds.size() == 1) { return triangleIds[0]; } - auto baseNorm = GetNormal(base); - size_t maxDotId = -1; double maxDot = -DBL_MAX; - for (auto id : triangleIds) { if (id == base.id) continue; - auto triNorm = GetNormal(A.triangles[id]); - auto other = A.triangles[id].GetNotShared(base); - auto above = CalcTriPt(base, other); - auto dot = (glm::dot(baseNorm, triNorm) + 1.0) / 2.0; // range dot from 0 to 1, positive for above, negative for below - if (above == TriangleVsPoint::BELOW) { dot = dot - 1.0; @@ -716,25 +586,20 @@ namespace fuzzybools { dot = 1.0 - dot; } - if (dot > maxDot) { maxDot = dot; maxDotId = id; } } - return maxDotId; } - //============================================================================================ - size_t GetPointWithMaxY() { double max = -DBL_MAX; size_t pointID = 0; - - for (auto &p : points) + for (auto& p : points) { if (p.location3D.y > max) { @@ -742,44 +607,34 @@ namespace fuzzybools pointID = p.id; } } - return pointID; } - //============================================================================================ - enum class TriangleVsPoint { ABOVE, BELOW, ON }; - //============================================================================================ - - TriangleVsPoint CalcTriPt(Triangle &T, size_t point) + TriangleVsPoint CalcTriPt(Triangle& T, size_t point) { auto norm = GetNormal(T); - auto dpt3d = points[point].location3D - points[T.a].location3D; auto dot = glm::dot(norm, dpt3d); - if (std::fabs(dot) < EPS_BIG) return TriangleVsPoint::ON; if (dot > 0.0) return TriangleVsPoint::ABOVE; return TriangleVsPoint::BELOW; } - //============================================================================================ - // simplify, see FindUppermostTriangleId - bool ShouldFlip(Triangle &T, Triangle &neighbour) + bool ShouldFlip(Triangle& T, Triangle& neighbour) { /* for each n in N, orient n by formula: if triangle T and N are vertices ABCDEF with BCDE as the shared edge, consider three cases: - E is above T, then A is above n if dot(normal(T), E) > 0 => dot(normal(n), A) > 0, else flip n E is below T, then A is below n @@ -789,19 +644,15 @@ namespace fuzzybools */ auto normT = GetNormal(T); auto normNB = GetNormal(neighbour); - if (T.SamePoints(neighbour)) { // same tri, different winding, flip if not the same normal return glm::dot(normT, normNB) < 1.0 - EPS_BIG; } - auto A = T.GetNotShared(neighbour); auto E = neighbour.GetNotShared(T); - TriangleVsPoint EvsT = CalcTriPt(T, E); - TriangleVsPoint AvsN = CalcTriPt(neighbour, E); - + TriangleVsPoint AvsN = CalcTriPt(neighbour, A); if (EvsT == TriangleVsPoint::ABOVE && AvsN != TriangleVsPoint::ABOVE) { return true; @@ -820,95 +671,95 @@ namespace fuzzybools return false; } } - //============================================================================================ - - Vec3 GetNormal(Triangle &tri) + Vec3 GetNormal(Triangle& tri) { Vec3 temp(-1.0, -1.0, -1.0); Vec3 norm = glm::normalize(temp); computeSafeNormal(points[tri.a].location3D, points[tri.b].location3D, points[tri.c].location3D, norm, EPS_SMALL); return norm; } - //============================================================================================ - - SegmentSet &GetSegSetA() + SegmentSet& GetSegSetA() { return A; } - //============================================================================================ - - size_t AddPoint(const Vec3 &newPoint) - { - for (auto &pt : points) - { - if (pt == newPoint) - { - return pt.id; + std::unordered_map> pointHashMap; // Maps quantized key to list of point IDs + // Optimized AddPoint: + size_t AddPoint(const Vec3& newPoint) { + // Quantize the point based on tolerance + double invTol = 1.0 / toleranceVectorEquality; + glm::ivec3 primaryKey( + static_cast(std::floor(newPoint.x * invTol)), + static_cast(std::floor(newPoint.y * invTol)), + static_cast(std::floor(newPoint.z * invTol)) + ); + // Check primary and neighboring cells with radius 2 + const int RADIUS = 2; + for (int dx = -RADIUS; dx <= RADIUS; ++dx) { + for (int dy = -RADIUS; dy <= RADIUS; ++dy) { + for (int dz = -RADIUS; dz <= RADIUS; ++dz) { + glm::ivec3 neighKey = primaryKey + glm::ivec3(dx, dy, dz); + auto it = pointHashMap.find(neighKey); + if (it != pointHashMap.end()) { + for (size_t existingId : it->second) { + if (equals(points[existingId].location3D, newPoint, toleranceVectorEquality)) { + return existingId; // Match found + } + } + } + } } } - + // No match: add new point Point p; p.id = points.size(); p.location3D = newPoint; - points.push_back(p); - + // Insert into primary cell's list + pointHashMap[primaryKey].push_back(p.id); return p.id; } - + //============================================================================================ - - size_t AddPlane(const Vec3 &normal, double d, uint32_t refId) + size_t AddPlane(const Vec3& normal, double d, uint32_t refId) { - for (auto &plane : planes) + for (auto& plane : planes) { if (plane.refPlane == refId || plane.IsEqualTo(normal, d)) { return plane.id; } } - Plane p; p.id = planes.size(); p.refPlane = refId; p.normal = glm::normalize(normal); p.distance = d; planes.push_back(p); - return p.id; } - //============================================================================================ - - void Construct(const Geometry &A, const Geometry &B, bool isUnion) + void Construct(const Geometry& A, const Geometry& B, bool isUnion) { auto boxA = A.GetAABB(); auto boxB = B.GetAABB(); - AddGeometry(A, B, boxB, true, isUnion, 0); AddGeometry(B, A, boxA, false, isUnion, A.planes.size()); - _linkedA = &A; _linkedB = &B; } - //============================================================================================ - - void AddGeometry(const Geometry &geom, const Geometry &secondGeom, const AABB &relevantBounds, bool isA, bool isUnion, uint32_t offsetPlane) + void AddGeometry(const Geometry& geom, const Geometry& secondGeom, const AABB& relevantBounds, bool isA, bool isUnion, uint32_t offsetPlane) { #ifdef CSG_DEBUG_OUTPUT Geometry relevant; #endif - for (size_t i = 0; i < geom.numFaces; i++) { Face f = geom.GetFace(i); - auto faceBox = geom.GetFaceBox(i); - if (!faceBox.intersects(relevantBounds)) { if (isA) @@ -919,23 +770,18 @@ namespace fuzzybools { B.irrelevantFaces.push_back(i); } - continue; } - bool contact = false; - for (size_t j = 0; j < secondGeom.numFaces; j++) { auto faceBox2 = secondGeom.GetFaceBox(j); - if (faceBox.intersects(faceBox2)) { contact = true; break; } } - if (!contact) { if (isA) @@ -949,38 +795,31 @@ namespace fuzzybools B.irrelevantFaces_toTest.push_back(i); } } - continue; } - if (isA) { #ifdef CSG_DEBUG_OUTPUT -// DumpGeometry(geom, L"Initial_A.obj"); + // DumpGeometry(geom, L"Initial_A.obj"); #endif } else { #ifdef CSG_DEBUG_OUTPUT -// DumpGeometry(geom, L"Initial_B.obj"); + // DumpGeometry(geom, L"Initial_B.obj"); #endif } - auto a = geom.GetPoint(f.i0); auto b = geom.GetPoint(f.i1); auto c = geom.GetPoint(f.i2); - #ifdef CSG_DEBUG_OUTPUT -// relevant.AddFace(a, b, c, -1); + // relevant.AddFace(a, b, c, -1); #endif - Vec3 norm; if (computeSafeNormal(a, b, c, norm, EPS_SMALL)) { double rs = glm::dot(geom.planes[f.pId].normal, norm); - size_t planeId = -1; - if (rs < 0) { planeId = AddPlane(-geom.planes[f.pId].normal, -geom.planes[f.pId].distance, f.pId + offsetPlane); @@ -989,15 +828,12 @@ namespace fuzzybools { planeId = AddPlane(geom.planes[f.pId].normal, geom.planes[f.pId].distance, f.pId + offsetPlane); } - auto ia = AddPoint(a); auto ib = AddPoint(b); auto ic = AddPoint(c); - double da = glm::dot(norm, a); double db = glm::dot(norm, b); double dc = glm::dot(norm, c); - if (!planes[planeId].IsPointOnPlane(a)) { if (messages) @@ -1034,7 +870,6 @@ namespace fuzzybools planes[planeId].AddPoint(a); planes[planeId].AddPoint(b); planes[planeId].AddPoint(c); - if (isA) { A.AddFace(planeId, ia, ib, ic); @@ -1054,7 +889,6 @@ namespace fuzzybools } } } - if (isA) { relevantBVHA = MakeBVH(relevantA); @@ -1063,49 +897,40 @@ namespace fuzzybools { relevantBVHB = MakeBVH(relevantB); } - #ifdef CSG_DEBUG_OUTPUT // if (isA) // { - // DumpGeometry(relevant, L"relevantA.obj"); + // DumpGeometry(relevant, L"relevantA.obj"); // } // else // { - // DumpGeometry(relevant, L"relevantB.obj"); + // DumpGeometry(relevant, L"relevantB.obj"); // } #endif } - //============================================================================================ - - std::vector GetPointsOnPlane(Plane &p) + std::vector GetPointsOnPlane(Plane& p) { auto cp = planeToPoints[p.id]; std::sort(cp.begin(), cp.end()); cp.erase(std::unique(cp.begin(), cp.end()), cp.end()); return cp; } - //============================================================================================ - // pair of lineID, distance - std::vector> BuildSegments(const std::vector &a, const std::vector &b) const + std::vector> BuildSegments(const std::vector& a, const std::vector& b) const { if (a.size() == 0 || b.size() == 0) { return {}; } - // we need to figure out the overlap between the two lists of intersections // we can be clever here and try to conclude that the first point must be the start of a segment and the next point would end that segment // however, we would really shoot ourselves in the foot as any coplanar results are missing from these two sets // let's just make some segments that span both intersection lists, and eat the overhead - double min = std::max(a[0], b[0]); double max = std::min(a[a.size() - 1], b[b.size() - 1]); - std::vector points; - for (size_t i = 0; i < a.size(); i++) { double val = a[i]; @@ -1114,7 +939,6 @@ namespace fuzzybools points.push_back(val); } } - for (size_t i = 0; i < b.size(); i++) { double val = b[i]; @@ -1123,92 +947,35 @@ namespace fuzzybools points.push_back(val); } } - const auto double_less = +[](double left, double right) - { return left < right; }; + { return left < right; }; std::sort(points.begin(), points.end(), double_less); - std::vector> result; result.reserve(points.size() * 2); - for (size_t i = 1; i < points.size(); i++) { result.emplace_back(points[i - 1], points[i]); } - // Remove redundant lines result.erase(std::unique(result.begin(), result.end()), result.end()); - return result; } - //============================================================================================ - - std::vector> GetNonIntersectingSegments(Line &l) + std::vector> GetNonIntersectingSegments(Line& l) { - std::vector> pointsInOrder; - - for (auto segment : l.GetSegments()) - { - if (!l.IsPointOnLine(points[segment.first].location3D)) - { - if (messages) - { - printf("point not on line in GetNonIntersectingSegments\n"); - } - if (messages) - { - printf("points[segment.first].location3D = (%12.8f, %12.8f, %12.8f)\n", points[segment.first].location3D.x, points[segment.first].location3D.y, points[segment.first].location3D.z); - } - if (messages) - { - printf("l.origin = (%12.8f, %12.8f, %12.8f)\n", l.origin.x, l.origin.y, l.origin.z); - } - if (messages) - { - printf("l.direction = (%12.8f, %12.8f, %12.8f)\n", l.direction.x, l.direction.y, l.direction.z); - } - } - - if (!l.IsPointOnLine(points[segment.second].location3D)) - { - if (messages) - { - printf("point not on line in GetNonIntersectingSegments\n"); - } - if (messages) - { - printf("points[segment.second].location3D = (%12.8f, %12.8f, %12.8f)\n", points[segment.second].location3D.x, points[segment.second].location3D.y, points[segment.second].location3D.z); - } - if (messages) - { - printf("l.origin = (%12.8f, %12.8f, %12.8f)\n", l.origin.x, l.origin.y, l.origin.z); - } - if (messages) - { - printf("l.direction = (%12.8f, %12.8f, %12.8f)\n", l.direction.x, l.direction.y, l.direction.z); - } - } - - pointsInOrder.emplace_back(segment.first, l.GetPosOnLine(points[segment.first].location3D)); - pointsInOrder.emplace_back(segment.second, l.GetPosOnLine(points[segment.second].location3D)); - } - + std::vector> pointsInOrder = l.points; std::sort( - pointsInOrder.begin(), pointsInOrder.end(), [&](const std::pair &left, const std::pair &right) - { return left.second > right.second; }); - + pointsInOrder.begin(), pointsInOrder.end(), [&](const std::pair& left, const std::pair& right) + { return left.first < right.first; }); std::vector> segmentsWithoutIntersections; - if (pointsInOrder.empty()) return {}; - - size_t cur = pointsInOrder[0].first; + size_t cur = pointsInOrder[0].second; for (size_t i = 1; i < pointsInOrder.size(); i++) { // this fills gaps, but we don't care about that in this stage // gaps filled will be removed during inside/outside checking - size_t next = pointsInOrder[i].first; + size_t next = pointsInOrder[i].second; if (cur != next) { segmentsWithoutIntersections.emplace_back(cur, next); @@ -1217,47 +984,35 @@ namespace fuzzybools } return segmentsWithoutIntersections; } - //============================================================================================ - - void TriangulatePlane(Geometry &geom, Plane &p) + void TriangulatePlane(Geometry& geom, Plane& p) { - // grab all points on the plane auto pointsOnPlane = GetPointsOnPlane(p); - // temporarily project all points for triangulation // NOTE: these points should not be used as output, only as placeholder for triangulation! - auto basis = p.MakeBasis(); - std::unordered_map pointToProjectedPoint; pointToProjectedPoint.reserve(pointsOnPlane.size()); - std::unordered_map projectedPointToPoint; projectedPointToPoint.reserve(pointsOnPlane.size()); - std::vector projectedPoints; projectedPoints.reserve(pointsOnPlane.size()); - - for (auto &pointId : pointsOnPlane) + for (auto& pointId : pointsOnPlane) { pointToProjectedPoint[pointId] = projectedPoints.size(); projectedPointToPoint[projectedPoints.size()] = pointId; projectedPoints.push_back(basis.project(points[pointId].location3D)); } - std::set> edges; std::set> defaultEdges; static int i = 0; i++; - - for (auto &line : p.lines) + for (auto& line : p.lines) { // these segments might intersect internally, lets resolve that so we get a valid chain auto segments = GetNonIntersectingSegments(line); - - for (auto &segment : segments) + for (auto& segment : segments) { if (pointToProjectedPoint.count(segment.first) == 0) { @@ -1266,7 +1021,6 @@ namespace fuzzybools { printf("unknown point in list, repairing in TriangulateLine\n"); } - pointToProjectedPoint[segment.first] = projectedPoints.size(); projectedPointToPoint[projectedPoints.size()] = segment.first; projectedPoints.push_back(basis.project(points[segment.first].location3D)); @@ -1278,220 +1032,164 @@ namespace fuzzybools { printf("unknown point in list, repairing in TriangulateLine\n"); } - pointToProjectedPoint[segment.second] = projectedPoints.size(); projectedPointToPoint[projectedPoints.size()] = segment.second; projectedPoints.push_back(basis.project(points[segment.second].location3D)); } - auto projectedIndexA = pointToProjectedPoint[segment.first]; auto projectedIndexB = pointToProjectedPoint[segment.second]; - if (projectedIndexA != projectedIndexB) { defaultEdges.insert(segment); - edges.insert(std::make_pair(projectedIndexA, projectedIndexB)); + edges.insert(std::make_pair(std::min(projectedIndexA, projectedIndexB), std::max(projectedIndexA, projectedIndexB))); } } } - #ifdef CSG_DEBUG_OUTPUT // std::vector> edgesPrinted; - // for (auto &e : edges) // { - // edgesPrinted.push_back({projectedPoints[e.first], projectedPoints[e.second]}); + // edgesPrinted.push_back({projectedPoints[e.first], projectedPoints[e.second]}); // } - // DumpSVGLines(edgesPrinted, L"poly.html"); #endif - CDT::Triangulation cdt(CDT::VertexInsertionOrder::AsProvided); std::vector cdt_edges; cdt_edges.reserve(edges.size()); - std::vector> cdt_verts; cdt_verts.reserve(projectedPoints.size()); - - for (auto &point : projectedPoints) + for (auto& point : projectedPoints) { cdt_verts.emplace_back(CDT::V2d::make(point.x, point.y)); } - - for (auto &edge : edges) + for (auto& edge : edges) { cdt_edges.emplace_back((uint32_t)edge.first, (uint32_t)edge.second); } - auto mapping = CDT::RemoveDuplicatesAndRemapEdges(cdt_verts, cdt_edges).mapping; - cdt.insertVertices(cdt_verts); cdt.insertEdges(cdt_edges); - cdt.eraseSuperTriangle(); - auto triangles = cdt.triangles; - // auto contourLoop = FindLargestEdgeLoop(projectedPoints, edges); - #ifdef CSG_DEBUG_OUTPUT // std::vector> edges3DTriangles; // std::set> edgesTriangles; // std::set> finalEdgesTriangles; #endif - - for (auto &tri : triangles) + for (auto& tri : triangles) { #ifdef CSG_DEBUG_OUTPUT // edgesTriangles.insert(std::make_pair(tri.vertices[0], tri.vertices[1])); // edgesTriangles.insert(std::make_pair(tri.vertices[1], tri.vertices[2])); // edgesTriangles.insert(std::make_pair(tri.vertices[0], tri.vertices[2])); #endif - size_t pointIdA = projectedPointToPoint[mapping[tri.vertices[0]]]; size_t pointIdB = projectedPointToPoint[mapping[tri.vertices[1]]]; size_t pointIdC = projectedPointToPoint[mapping[tri.vertices[2]]]; - auto ptA = points[pointIdA].location3D; auto ptB = points[pointIdB].location3D; auto ptC = points[pointIdC].location3D; - glm::dvec3 v1 = glm::normalize(ptA - ptB); glm::dvec3 v2 = glm::normalize(ptA - ptC); glm::dvec3 v3 = glm::normalize(ptB - ptC); double rs1 = glm::dot(v1, v2); double rs2 = glm::dot(v2, v3); double rs3 = glm::dot(v1, v3); - if (std::abs(rs1) > 1 - toleranceThinTriangle || std::abs(rs2) > 1 - toleranceThinTriangle || std::abs(rs3) > 1 - toleranceThinTriangle) { continue; } - auto pt2DA = projectedPoints[mapping[tri.vertices[0]]]; auto pt2DB = projectedPoints[mapping[tri.vertices[1]]]; auto pt2DC = projectedPoints[mapping[tri.vertices[2]]]; - auto triCenter = (ptA + ptB + ptC) / 3.0; - Vec raydir = computeNormal(ptA, ptB, ptC); - auto posA = isInsideMesh(triCenter, glm::dvec3(0), relevantA, relevantBVHA, raydir); auto posB = isInsideMesh(triCenter, glm::dvec3(0), relevantB, relevantBVHB, raydir); - if (posA.loc != MeshLocation::BOUNDARY && posB.loc != MeshLocation::BOUNDARY) { continue; } - // If the 2D triangle is not inside the boundaries of the projected boundary of the face it requires further verification // It can't be discarded because inside/outside could fail when boundaries have internal partitions // Therefore new tests are required to verify that the triangle is on the boundary of A or B - glm::dvec2 t1 = projectedPoints[tri.vertices[0]]; glm::dvec2 t2 = projectedPoints[tri.vertices[1]]; glm::dvec2 t3 = projectedPoints[tri.vertices[2]]; - bool inside2d = isInsideBoundary(t1, t2, t3, edges, projectedPoints); - if (!inside2d) { - auto postA = isInsideMesh(triCenter, glm::dvec3(0), relevantA, relevantBVHA, raydir); auto postB = isInsideMesh(triCenter, glm::dvec3(0), relevantB, relevantBVHB, raydir); - if (postA.loc != MeshLocation::BOUNDARY && postB.loc != MeshLocation::BOUNDARY) { continue; } - auto ptt = glm::mix(triCenter, ptA, triangleEvaluationFactor); - postA = isInsideMesh(ptt, glm::dvec3(0), relevantA, relevantBVHA, raydir); postB = isInsideMesh(ptt, glm::dvec3(0), relevantB, relevantBVHB, raydir); - if (postA.loc != MeshLocation::BOUNDARY && postB.loc != MeshLocation::BOUNDARY) { continue; } - ptt = glm::mix(triCenter, ptB, triangleEvaluationFactor); - postA = isInsideMesh(ptt, glm::dvec3(0), relevantA, relevantBVHA, raydir); postB = isInsideMesh(ptt, glm::dvec3(0), relevantB, relevantBVHB, raydir); - if (postA.loc != MeshLocation::BOUNDARY && postB.loc != MeshLocation::BOUNDARY) { continue; } - ptt = glm::mix(triCenter, ptC, triangleEvaluationFactor); - postA = isInsideMesh(ptt, glm::dvec3(0), relevantA, relevantBVHA, raydir); postB = isInsideMesh(ptt, glm::dvec3(0), relevantB, relevantBVHB, raydir); - if (postA.loc != MeshLocation::BOUNDARY && postB.loc != MeshLocation::BOUNDARY) { continue; } } - // although CDT is great, it spits out too many or too little tris, we fix it manually // if (!IsPointInsideLoop(projectedPoints, contourLoop, triCenter)) //{ - // printf("removing point outside loop\n"); - // continue; + // printf("removing point outside loop\n"); + // continue; //} - // TODO: why is this swapped? winding doesnt matter much, but still geom.AddFace(ptB, ptA, ptC, p.refPlane); - #ifdef CSG_DEBUG_OUTPUT // edges3DTriangles.push_back({ glm::dvec2(ptA.z+ ptA.x/2, ptA.y+ ptA.x/2), glm::dvec2(ptB.z+ ptB.x/2, ptB.y+ ptB.x/2) }); // edges3DTriangles.push_back({ glm::dvec2(ptA.z+ ptA.x/2, ptA.y+ ptA.x/2), glm::dvec2(ptC.z+ ptC.x/2, ptC.y+ ptC.x/2) }); // edges3DTriangles.push_back({ glm::dvec2(ptB.z+ ptB.x/2, ptB.y+ ptB.x/2), glm::dvec2(ptC.z+ ptC.x/2, ptC.y+ ptC.x/2) }); // DumpSVGLines(edges3DTriangles, L"edges_tri.html"); #endif - #ifdef CSG_DEBUG_OUTPUT // finalEdgesTriangles.insert(std::make_pair(tri.vertices[0], tri.vertices[1])); // finalEdgesTriangles.insert(std::make_pair(tri.vertices[1], tri.vertices[2])); // finalEdgesTriangles.insert(std::make_pair(tri.vertices[0], tri.vertices[2])); #endif } - #ifdef CSG_DEBUG_OUTPUT // std::vector> edgesPrinted2; - // for (auto& e : edgesTriangles) // { - // edgesPrinted2.push_back({ projectedPoints[e.first], projectedPoints[e.second] }); + // edgesPrinted2.push_back({ projectedPoints[e.first], projectedPoints[e.second] }); // } - // DumpSVGLines(edgesPrinted2, L"poly_triangulation.html"); - // std::vector> finalEdgesPrinted; - // for (auto& e : finalEdgesTriangles) // { - // finalEdgesPrinted.push_back({ projectedPoints[e.first], projectedPoints[e.second] }); + // finalEdgesPrinted.push_back({ projectedPoints[e.first], projectedPoints[e.second] }); // } - // DumpSVGLines(finalEdgesPrinted, L"final_poly_triangulation.html"); #endif } //============================================================================================ - std::unordered_map> planeToLines; - //============================================================================================ - std::unordered_map> planeToPoints; - //============================================================================================ - void AddRefPlaneToPoint(size_t point, size_t plane) { ReferencePlane ref; @@ -1501,19 +1199,15 @@ namespace fuzzybools planeToPoints[plane].push_back(point); } }; - //============================================================================================ - - inline void AddSegments(Plane &p, SharedPosition &sp, Line &templine, const std::vector> &segments) + inline void AddSegments(Plane& p, SharedPosition& sp, Line& templine, const std::vector>& segments) { // NOTE: this is a design flaw, the addline may return a line that is // EQUIVALENT BUT NOT IDENTICAL // hence line distances mentioned in "segments" DO apply to templine // but not necessary (but possibly) to isectLineId auto isectLineId = p.AddLine(templine.origin, templine.direction); - - auto &isectLine = p.lines[isectLineId.first]; - + auto& isectLine = p.lines[isectLineId.first]; if (!p.IsPointOnPlane(isectLine.origin) || !p.IsPointOnPlane(isectLine.origin + isectLine.direction * 100.)) { if (messages) @@ -1521,11 +1215,9 @@ namespace fuzzybools printf("Bad isect line in AddSegments\n"); } } - - for (auto &seg : segments) + for (auto& seg : segments) { auto pos = templine.GetPosOnLine(seg.first); - if (!p.aabb.contains(pos)) { if (messages) @@ -1533,10 +1225,8 @@ namespace fuzzybools printf("making pos outside in AddSegments]\n"); } } - size_t ptA = sp.AddPoint(pos); size_t ptB = sp.AddPoint(templine.GetPosOnLine(seg.second)); - if (!p.aabb.contains(sp.points[ptA].location3D)) { if (messages) @@ -1551,13 +1241,10 @@ namespace fuzzybools printf("bad points in AddS.segments\n"); } } - - // if (ptA != ptB) { isectLine.AddPointToLine(isectLine.GetPosOnLine(sp.points[ptA].location3D), ptA); isectLine.AddPointToLine(isectLine.GetPosOnLine(sp.points[ptB].location3D), ptB); } - if (!p.IsPointOnPlane(sp.points[ptA].location3D)) { if (messages) @@ -1572,167 +1259,84 @@ namespace fuzzybools printf("bad point in AddSegments\n"); } } - sp.AddRefPlaneToPoint(ptA, p.id); sp.AddRefPlaneToPoint(ptB, p.id); } } - //============================================================================================ - - inline std::vector ComputeInitialIntersections(Plane &p, SharedPosition &sp, const Line &lineA) + inline std::vector ComputeInitialIntersections(Plane& p, SharedPosition& sp, const Line& lineA) { - double size = 1.0E+08; // TODO: this is bad - - for (auto &point : sp.points) - { - const auto d2 = glm::distance2(lineA.origin, point.location3D); - size = std::max(size, d2); - } - - size = std::sqrt(size); - - auto Astart = lineA.origin + lineA.direction * (size * 2); - auto Aend = lineA.origin - lineA.direction * (size * 2); - std::vector distances; distances.reserve(p.lines.size()); - // line B is expected to have the segments already filled, line A is not - for (auto &line : p.lines) + for (auto& line : p.lines) { // skip collinear if (lineA.IsCollinear(line)) continue; - - for (const auto &seg : line.GetSegments()) + Vec3 relative = lineA.origin - line.origin; + Vec3 cr = glm::cross(lineA.direction, line.direction); + double denom = glm::dot(cr, p.normal); + if (std::fabs(denom) < EPS_SMALL) + continue; + double t = -glm::dot(glm::cross(relative, line.direction), p.normal) / denom; + Vec3 point = lineA.origin + t * lineA.direction; + if (!line.IsPointOnLine(point)) + continue; + double distB = line.GetPosOnLine(point); + if (!line.ContainsPos(distB)) + continue; + if (!p.IsPointOnPlane(point)) { - auto result = LineLineIntersection( - Astart, - Aend, - sp.points[seg.first].location3D, - sp.points[seg.second].location3D); - - if (result.distance < _TOLERANCE_PLANE_INTERSECTION) + if (messages) { - if (!p.aabb.contains(sp.points[seg.first].location3D)) - { - if (messages) - { - printf("bad points in ComputeInitialIntersections\n"); - } - } - if (!p.aabb.contains(sp.points[seg.second].location3D)) - { - if (messages) - { - printf("bad points in ComputeInitialIntersections\n"); - } - } - - // intersection, mark index of line B and distance on line A - distances.emplace_back(lineA.GetPosOnLine(result.point2)); - auto pt = lineA.GetPosOnLine(distances[distances.size() - 1]); - - if (!p.aabb.contains(result.point2)) - { - if (messages) - { - printf("bad points in ComputeInitialIntersections\n"); - } - } - - if (!equals(pt, result.point2, _TOLERANCE_PLANE_INTERSECTION)) - { - if (messages) - { - printf("BAD POINT in ComputeInitialIntersections\n"); - } - } + printf("bad point in ComputeInitialIntersections\n"); } } + double distA = lineA.GetPosOnLine(point); + distances.emplace_back(distA); } - const auto double_less = +[](double left, double right) - { return left < right; }; + { return left < right; }; std::sort(distances.begin(), distances.end(), double_less); - distances.erase(std::unique(distances.begin(), distances.end()), distances.end()); - return distances; } - //============================================================================================ - - inline void AddLineLineIntersections(Plane &p, SharedPosition &sp, Line &lineA, Line &lineB) + inline void AddLineLineIntersections(Plane& p, SharedPosition& sp, Line& lineA, Line& lineB) { - for (auto segA : lineA.GetSegments()) - { - for (auto segB : lineB.GetSegments()) - { - // check isect A vs B - if (!p.HasOverlap(segA, segB)) - { - // no overlap, possibility of intersection - auto result = LineLineIntersection( - sp.points[segA.first].location3D, - sp.points[segA.second].location3D, - sp.points[segB.first].location3D, - sp.points[segB.second].location3D); - - if (result.distance < SCALED_EPS_BIG) - { - // intersection! Take center and insert - if (!p.aabb.contains(result.point1)) - { - if (messages) - { - printf("bad points in AddLineLineIntersections\n"); - } - continue; - } - - size_t point = sp.AddPoint((result.point1)); - - lineA.AddPointToLine(lineA.GetPosOnLine(sp.points[point].location3D), point); - lineB.AddPointToLine(lineB.GetPosOnLine(sp.points[point].location3D), point); - - // point falls on intersection of two lines, so is part of those lines and all planes of those lines - { - ReferenceLine ref; - ref.pointID = point; - ref.lineID = lineA.id; - ref.location = lineA.GetPosOnLine(result.point1); - sp.points[point].lines.push_back(ref); - - // TODO: FIX THIS POINT MIGHT NOT BE ON THE PLANE ACTUALLY - - for (auto &plane : lineA.planes) - { - sp.AddRefPlaneToPoint(point, plane.planeID); - } - } - { - ReferenceLine ref; - ref.pointID = point; - ref.lineID = lineB.id; - ref.location = lineB.GetPosOnLine(result.point2); - sp.points[point].lines.push_back(ref); - - for (auto &plane : lineB.planes) - { - sp.AddRefPlaneToPoint(point, plane.planeID); - } - } - } - } - } - } + if (lineA.IsCollinear(lineB)) + return; + Vec3 relative = lineA.origin - lineB.origin; + Vec3 cr = glm::cross(lineA.direction, lineB.direction); + double denom = glm::dot(cr, p.normal); + if (std::fabs(denom) < EPS_SMALL) + return; + double t = -glm::dot(glm::cross(relative, lineB.direction), p.normal) / denom; + Vec3 point = lineA.origin + t * lineA.direction; + if (!lineA.IsPointOnLine(point) || !lineB.IsPointOnLine(point)) + return; + double distA = lineA.GetPosOnLine(point); + double distB = lineB.GetPosOnLine(point); + if (!lineA.ContainsPos(distA) || !lineB.ContainsPos(distB)) + return; + size_t point_id = sp.AddPoint(point); + lineA.AddPointToLine(distA, point_id); + lineB.AddPointToLine(distB, point_id); + ReferenceLine refA; + refA.pointID = point_id; + refA.lineID = lineA.id; + refA.location = distA; + sp.points[point_id].lines.push_back(refA); + ReferenceLine refB; + refB.pointID = point_id; + refB.lineID = lineB.id; + refB.location = distB; + sp.points[point_id].lines.push_back(refB); + sp.AddRefPlaneToPoint(point_id, p.id); } - //============================================================================================ - - inline void AddLineLineIsects(Plane &p, SharedPosition &sp) + inline void AddLineLineIsects(Plane& p, SharedPosition& sp) { for (size_t lineAIndex = 0; lineAIndex < p.lines.size(); lineAIndex++) { @@ -1742,65 +1346,50 @@ namespace fuzzybools } } } - //============================================================================================ - - inline Geometry Normalize(const Geometry &A, const Geometry &B, SharedPosition &sp, bool UNION) + inline Geometry Normalize(const Geometry& A, const Geometry& B, SharedPosition& sp, bool UNION) { - // construct all contours, derive lines auto contoursA = sp.A.GetContourSegments(); - - for (auto &[planeId, contours] : contoursA) + for (auto& [planeId, contours] : contoursA) { std::vector> edges; - - Plane &p = sp.planes[planeId]; - + Plane& p = sp.planes[planeId]; #ifdef CSG_DEBUG_OUTPUT // auto basis = p.MakeBasis(); - // for (auto& segment : contours) // { - // edges.push_back({ basis.project(sp.points[segment.first].location3D), basis.project(sp.points[segment.second].location3D) }); + // edges.push_back({ basis.project(sp.points[segment.first].location3D), basis.project(sp.points[segment.second].location3D) }); // } // DumpSVGLines(edges, L"contour_A.html"); #endif - - for (auto &segment : contours) + for (auto& segment : contours) { auto lineId = sp.planes[planeId].AddLine(sp.points[segment.first], sp.points[segment.second]); } } - auto contoursB = sp.B.GetContourSegments(); - - for (auto &[planeId, contours] : contoursB) + for (auto& [planeId, contours] : contoursB) { std::vector> edges; - - Plane &p = sp.planes[planeId]; - + Plane& p = sp.planes[planeId]; #ifdef CSG_DEBUG_OUTPUT // auto basis = p.MakeBasis(); - // for (auto& segment : contours) // { - // edges.push_back({ basis.project(sp.points[segment.first].location3D), basis.project(sp.points[segment.second].location3D) }); + // edges.push_back({ basis.project(sp.points[segment.first].location3D), basis.project(sp.points[segment.second].location3D) }); // } // DumpSVGLines(edges, L"contour_B.html"); #endif - - for (auto &segment : contours) + for (auto& segment : contours) { auto lineId = sp.planes[planeId].AddLine(sp.points[segment.first], sp.points[segment.second]); } } - // put all points on lines/planes - for (auto &p : sp.points) + for (auto& p : sp.points) { - for (auto &plane : sp.planes) + for (auto& plane : sp.planes) { if (plane.IsPointOnPlane(p.location3D)) { @@ -1809,43 +1398,31 @@ namespace fuzzybools } } } - - for (auto &plane : sp.planes) - { - AddLineLineIsects(plane, sp); - } - // intersect planes for (size_t planeAIndex = 0; planeAIndex < sp.planes.size(); planeAIndex++) { - for (size_t planeBIndex = 0; planeBIndex < sp.planes.size(); planeBIndex++) + for (size_t planeBIndex = planeAIndex + 1; planeBIndex < sp.planes.size(); planeBIndex++) { - auto &planeA = sp.planes[planeAIndex]; - auto &planeB = sp.planes[planeBIndex]; - + auto& planeA = sp.planes[planeAIndex]; + auto& planeB = sp.planes[planeBIndex]; if (!planeA.aabb.intersects(planeB.aabb)) { continue; } - // plane intersect results in new lines // new lines result in new line intersects // new line intersects result in new points - if (std::fabs(glm::dot(planeA.normal, planeB.normal)) > 1.0 - EPS_BIG) { // parallel planes, don't care continue; } - // calculate plane intersection line auto result = PlanePlaneIsect(planeA.normal, planeA.distance, planeB.normal, planeB.distance); - // TODO: invalid temp line object Line intersectionLine; intersectionLine.origin = result.pos; intersectionLine.direction = result.dir; - if (!planeA.IsPointOnPlane(intersectionLine.origin) || !planeA.IsPointOnPlane(intersectionLine.origin + intersectionLine.direction * 1000.)) { if (messages) @@ -1860,14 +1437,11 @@ namespace fuzzybools printf("Bad isect line in Normalize\n"); } } - // get all intersection points with the shared line and both planes auto isectA = ComputeInitialIntersections(planeA, sp, intersectionLine); auto isectB = ComputeInitialIntersections(planeB, sp, intersectionLine); - // from these, figure out the shared segments on the current line produced by these two planes auto segments = sp.BuildSegments(isectA, isectB); - if (segments.empty()) { // nothing resulted from this plane-plane intersection @@ -1880,15 +1454,13 @@ namespace fuzzybools } } } - - for (auto &plane : sp.planes) + for (auto& plane : sp.planes) { AddLineLineIsects(plane, sp); } - - for (auto &p : sp.points) + for (auto& p : sp.points) { - for (auto &plane : sp.planes) + for (auto& plane : sp.planes) { if (plane.IsPointOnPlane(p.location3D)) { @@ -1896,53 +1468,42 @@ namespace fuzzybools } } } - // from the inserted geometries, all lines planes and points are now merged into a single set of shared planes lines and points // from this starting point, we can triangulate all planes and obtain the triangulation of the intersected set of geometries // this mesh itself is not a boolean result, but rather a merging of all operands - Geometry geom; - for (auto &plane : sp.planes) + for (auto& plane : sp.planes) { - #ifdef CSG_DEBUG_OUTPUT // std::vector> edges; - // auto basis = plane.MakeBasis(); - // for (auto& line : plane.lines) { - // // Get line parameters - // auto origin = line.origin; // 3D point (glm::dvec3) - // auto direction = line.direction; // 3D vector (glm::dvec3) - - // // Convert each distance to a 3D point along the line - // std::vector lineSegments; - // for (auto distance : line.points) { - // glm::dvec3 point3D = origin + glm::dvec3(direction.x * distance.first, direction.y * distance.first, direction.z * distance.first); // 3D calculation - // glm::dvec2 point2D = basis.project(point3D); // Project to 2D - // lineSegments.push_back(point2D); - // } - - // // Create edges between consecutive points - // for (size_t i = 0; i < lineSegments.size() - 1; ++i) { - // edges.push_back({ - // lineSegments[i], - // lineSegments[i + 1] - // }); - // } + // // Get line parameters + // auto origin = line.origin; // 3D point (glm::dvec3) + // auto direction = line.direction; // 3D vector (glm::dvec3) + // // Convert each distance to a 3D point along the line + // std::vector lineSegments; + // for (auto distance : line.points) { + // glm::dvec3 point3D = origin + glm::dvec3(direction.x * distance.first, direction.y * distance.first, direction.z * distance.first); // 3D calculation + // glm::dvec2 point2D = basis.project(point3D); // Project to 2D + // lineSegments.push_back(point2D); + // } + // // Create edges between consecutive points + // for (size_t i = 0; i < lineSegments.size() - 1; ++i) { + // edges.push_back({ + // lineSegments[i], + // lineSegments[i + 1] + // }); + // } // } - // DumpSVGLines(edges, L"contour.html"); #endif - sp.TriangulatePlane(geom, plane); - #ifdef CSG_DEBUG_OUTPUT // DumpGeometry(geom, L"triangulated.obj"); #endif } - - for (auto &plane : A.planes) + for (auto& plane : A.planes) { SimplePlane p; p.normal = plane.normal; @@ -1950,8 +1511,7 @@ namespace fuzzybools geom.planes.push_back(p); geom.hasPlanes = true; } - - for (auto &plane : B.planes) + for (auto& plane : B.planes) { SimplePlane p; p.normal = plane.normal; @@ -1959,60 +1519,45 @@ namespace fuzzybools geom.planes.push_back(p); geom.hasPlanes = true; } - uint32_t offsetA = A.planes.size(); - // re-add irrelevant faces that should be tested - for (auto &faceIndex : sp.A.irrelevantFaces_toTest) + for (auto& faceIndex : sp.A.irrelevantFaces_toTest) { - const Face &f = sp._linkedA->GetFace(faceIndex); - + const Face& f = sp._linkedA->GetFace(faceIndex); auto a = sp._linkedA->GetPoint(f.i0); auto b = sp._linkedA->GetPoint(f.i1); auto c = sp._linkedA->GetPoint(f.i2); - geom.AddFace(a, b, c, f.pId); } - - for (auto &faceIndex : sp.B.irrelevantFaces_toTest) + for (auto& faceIndex : sp.B.irrelevantFaces_toTest) { - const Face &f = sp._linkedB->GetFace(faceIndex); - + const Face& f = sp._linkedB->GetFace(faceIndex); auto a = sp._linkedB->GetPoint(f.i0); auto b = sp._linkedB->GetPoint(f.i1); auto c = sp._linkedB->GetPoint(f.i2); - geom.AddFace(a, b, c, f.pId + offsetA); } - geom.data = geom.numFaces; - // re-add irrelevant faces - for (auto &faceIndex : sp.A.irrelevantFaces) + for (auto& faceIndex : sp.A.irrelevantFaces) { - const Face &f = sp._linkedA->GetFace(faceIndex); - + const Face& f = sp._linkedA->GetFace(faceIndex); auto a = sp._linkedA->GetPoint(f.i0); auto b = sp._linkedA->GetPoint(f.i1); auto c = sp._linkedA->GetPoint(f.i2); - geom.AddFace(a, b, c, f.pId); } - if (UNION) { - for (auto &faceIndex : sp.B.irrelevantFaces) + for (auto& faceIndex : sp.B.irrelevantFaces) { - const Face &f = sp._linkedB->GetFace(faceIndex); - + const Face& f = sp._linkedB->GetFace(faceIndex); auto a = sp._linkedB->GetPoint(f.i0); auto b = sp._linkedB->GetPoint(f.i1); auto c = sp._linkedB->GetPoint(f.i2); - geom.AddFace(a, b, c, f.pId + offsetA); } } - return geom; } -} +} \ No newline at end of file diff --git a/src/cpp/web-ifc/geometry/operations/geometryutils.h b/src/cpp/web-ifc/geometry/operations/geometryutils.h index 8dd57475..14a2defa 100644 --- a/src/cpp/web-ifc/geometry/operations/geometryutils.h +++ b/src/cpp/web-ifc/geometry/operations/geometryutils.h @@ -16,6 +16,7 @@ #include #include #include +#include #include "../representation/geometry.h" #include "../representation/IfcGeometry.h" #include @@ -117,76 +118,61 @@ namespace webifc::geometry return bimGeometry::computeSafeNormal(v1, v2, v3, normal, eps); } - inline bool GetBasisFromCoplanarPoints(std::vector &points, glm::dvec3 &v1, glm::dvec3 &v2, glm::dvec3 &v3) + inline bool GetBasisFromCoplanarPoints( + const std::vector& points, + glm::dvec3& v1, + glm::dvec3& v2, + glm::dvec3& v3 + ) { + if (points.size() < 3) { + // Need at least three points to define a plane (basis) + return false; + } + v1 = points[0]; - for (auto &p : points) - { - if (v1 != p) - { - v2 = p; + // --- 1. Find the second distinct point (v2) --- + bool v2_found = false; + for (size_t i = 1; i < points.size(); ++i) { + if (glm::distance2(v1, points[i]) > EPS_SMALL * EPS_SMALL) { + v2 = points[i]; + v2_found = true; break; } } - - glm::dvec3 normal; - // multiple tries to find the best match - for (double i = 0; i < 4; i++) - { - double EPS = EPS_SMALL; - if (i == 0) - { - EPS = 100; - } - if (i == 1) - { - EPS = 1; - } - if (i == 2) - { - EPS = 0.01; - } - if (i == 3) - { - EPS = 1e-03; - } - for (auto &p : points) - { - if (computeSafeNormal(v1, v2, p, normal, EPS)) - { - v3 = p; - return true; - } - } + if (!v2_found) { + // All points are identical + return false; } - double d1 = 0; + // --- 2. Search for the third non-collinear point (v3) --- + glm::dvec3 v1v2 = v2 - v1; - for (auto &p : points) - { - double d2 = glm::distance(v1, p); - if (d1 < d2) - { - d1 = d2; - v2 = p; + // Use a robust check for collinearity: check the area of the triangle (v1, v2, p) + for (size_t i = 1; i < points.size(); ++i) { + const glm::dvec3& p = points[i]; + + // Skip v1 and v2, though the check below should handle them correctly + if (glm::distance2(v1, p) < EPS_SMALL * EPS_SMALL || glm::distance2(v2, p) < EPS_SMALL * EPS_SMALL) { + continue; } - } - d1 = 0; - for (auto &p : points) - { - double d2 = glm::distance(v1, p); - double d3 = glm::distance(v2, p); - if (d1 < d2 + d3) - { - d1 = d2 + d3; + glm::dvec3 v1p = p - v1; + + // Check if the cross product is significant (meaning the points are not collinear) + glm::dvec3 cross_product = glm::cross(v1v2, v1p); + + // Use a small, fixed tolerance (e.g., 1e-12) for the square of the cross product magnitude. + if (glm::length2(cross_product) > 1e-12) { v3 = p; + return true; } } - // multiple tries to find the best match - return computeSafeNormal(v1, v2, v3, normal, 1e-08); + // --- 3. Fallback: No three non-collinear points found --- + // All points are either identical or collinear. + return false; } inline void TriangulateBounds(IfcGeometry& geometry, std::vector& bounds, uint32_t expressID)