diff --git a/src/ifcgeom/IfcFace.cpp b/src/ifcgeom/IfcFace.cpp index 1c615d032d..fba37e4092 100644 --- a/src/ifcgeom/IfcFace.cpp +++ b/src/ifcgeom/IfcFace.cpp @@ -38,6 +38,7 @@ #include "../ifcgeom/IfcGeom.h" #include "../ifcgeom_schema_agnostic/face_definition.h" +#include "../ifcgeom_schema_agnostic/wire_utils.h" #define Kernel MAKE_TYPE_NAME(Kernel) @@ -156,7 +157,7 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcFace* l, TopoDS_Shape& result) } } else { gp_Pln pln; - if (approximate_plane_through_wire(wire, pln)) { + if (util::approximate_plane_through_wire(wire, pln, getValue(GV_PRECISION))) { fd.surface() = new Geom_Plane(pln); } } @@ -185,11 +186,17 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcFace* l, TopoDS_Shape& result) if (fd.all_outer()) { for (const auto& w : fd.wires()) { TopTools_ListOfShape fl; - triangulate_wire({ w }, fl); + auto r = util::triangulate_wire({ w }, fl); face_list.Append(fl); + if (faceset_helper_ && r == util::TRIANGULATE_WIRE_NON_MANIFOLD) { + faceset_helper_->non_manifold() = true; + } } } else { - triangulate_wire(fd.wires(), face_list); + auto r = util::triangulate_wire(fd.wires(), face_list); + if (faceset_helper_ && r == util::TRIANGULATE_WIRE_NON_MANIFOLD) { + faceset_helper_->non_manifold() = true; + } } } else if (!fd.all_outer()) { BRepBuilderAPI_MakeFace mf(fd.surface(), fd.outer_wire()); diff --git a/src/ifcgeom/IfcGeom.cpp b/src/ifcgeom/IfcGeom.cpp index 9e649bf015..d58a5aeacd 100644 --- a/src/ifcgeom/IfcGeom.cpp +++ b/src/ifcgeom/IfcGeom.cpp @@ -40,10 +40,6 @@ #include #include -#include -#include - - #include #include @@ -143,6 +139,7 @@ #include "../ifcgeom_schema_agnostic/IfcGeomTree.h" #include "../ifcgeom_schema_agnostic/boolean_utils.h" +#include "../ifcgeom_schema_agnostic/wire_utils.h" #include #include @@ -217,10 +214,6 @@ void MAKE_INIT_FN(KernelImplementation_)(IfcGeom::impl::KernelFactoryImplementat #define Kernel MAKE_TYPE_NAME(Kernel) -namespace { - -} - void IfcGeom::Kernel::set_offset(const std::array &p_offset) { offset = gp_Vec(p_offset[0], p_offset[1], p_offset[2]); @@ -897,13 +890,24 @@ bool IfcGeom::Kernel::convert_openings_fast(const IfcSchema::IfcProduct* entity, } #endif +double IfcGeom::Kernel::get_wire_intersection_tolerance(const TopoDS_Wire& wire) const { + return getValue(GV_NO_WIRE_INTERSECTION_TOLERANCE) > 0. + ? 0 + : faceset_helper_ + // eps is added to both ends of the parametric domain, so 3. is chosen to be on the safe side here. + ? (faceset_helper_->epsilon() / 3.) + // @todo re-evaluate 2. here for the reasons above: + : (std::min)(util::min_edge_length(wire) / 2., getValue(GV_PRECISION) * 10.); +} + bool IfcGeom::Kernel::convert_wire_to_face(const TopoDS_Wire& w, TopoDS_Face& face) { TopoDS_Wire wire = w; TopTools_ListOfShape results; - if (wire_intersections(wire, results)) { + + if (getValue(GV_NO_WIRE_INTERSECTION_CHECK) == 0. && util::wire_intersections(wire, results, get_wire_intersection_tolerance(wire), getValue(GV_PRECISION))) { Logger::Warning("Self-intersections with " + boost::lexical_cast(results.Extent()) + " cycles detected"); - select_largest(results, wire); + util::select_largest(results, wire); } bool is_2d = true; @@ -958,7 +962,7 @@ bool IfcGeom::Kernel::convert_wire_to_faces(const TopoDS_Wire& w, TopoDS_Compoun } TopTools_ListOfShape results; - if (wire_intersections(w, results)) { + if (getValue(GV_NO_WIRE_INTERSECTION_CHECK) == 0. && util::wire_intersections(w, results, get_wire_intersection_tolerance(w), getValue(GV_PRECISION))) { Logger::Warning("Self-intersections with " + boost::lexical_cast(results.Extent()) + " cycles detected"); } else { results.Clear(); @@ -3286,295 +3290,6 @@ bool IfcGeom::Kernel::is_identity_transform(IfcUtil::IfcBaseInterface* l) { } } -bool IfcGeom::Kernel::approximate_plane_through_wire(const TopoDS_Wire& wire, gp_Pln& plane, double eps) { - // Newell's Method is used for the normal calculation - // as a simple edge cross product can give opposite results - // for a concave face boundary. - // Reference: Graphics Gems III p. 231 - - const double eps_ = eps < 1. ? getValue(GV_PRECISION) : eps; - const double eps2 = eps_ * eps_; - - double x = 0, y = 0, z = 0; - gp_Pnt current, previous, first; - gp_XYZ center; - int n = 0; - - BRepTools_WireExplorer exp(wire); - - for (;; exp.Next()) { - const bool has_more = exp.More() != 0; - if (has_more) { - const TopoDS_Vertex& v = exp.CurrentVertex(); - current = BRep_Tool::Pnt(v); - center += current.XYZ(); - } else { - current = first; - } - if (n) { - const double& xn = previous.X(); - const double& yn = previous.Y(); - const double& zn = previous.Z(); - const double& xn1 = current.X(); - const double& yn1 = current.Y(); - const double& zn1 = current.Z(); - x += (yn - yn1)*(zn + zn1); - y += (xn + xn1)*(zn - zn1); - z += (xn - xn1)*(yn + yn1); - } else { - first = current; - } - if (!has_more) { - break; - } - previous = current; - ++n; - } - - if (n < 3) { - return false; - } - - plane = gp_Pln(center / n, gp_Dir(x, y, z)); - - exp.Init(wire); - for (; exp.More(); exp.Next()) { - const TopoDS_Vertex& v = exp.CurrentVertex(); - current = BRep_Tool::Pnt(v); - if (plane.SquareDistance(current) > eps2) { - return false; - } - } - - return true; -} - -bool IfcGeom::Kernel::flatten_wire(TopoDS_Wire& wire) { - gp_Pln pln; - if (!approximate_plane_through_wire(wire, pln)) { - return false; - } - TopoDS_Face face = BRepBuilderAPI_MakeFace(pln).Face(); - BRepAlgo_NormalProjection proj(face); - proj.Add(wire); - proj.Build(); - if (!proj.IsDone()) { - return false; - } - TopTools_ListOfShape list; - proj.BuildWire(list); - if (list.Extent() != 1) { - return false; - } - wire = TopoDS::Wire(list.First()); - return true; -} - -bool IfcGeom::Kernel::triangulate_wire(const std::vector& wires, TopTools_ListOfShape& faces) { - // This is a bit of a precarious approach, but seems to work for the - // versions of OCCT tested for. OCCT has a Delaunay triangulation function - // BRepMesh_Delaun, but it is notoriously hard to interpret the results - // (due to the Bowyer-Watson super triangle perhaps?). Therefore - // alternatively we use the regular OCCT incremental mesher on a new face - // created from the UV coordinates of the original wire. Pray to our gods - // that the vertex coordinates are unaffected by the meshing algorithm and - // map them back to 3d coordinates when iterating over the mesh triangles. - - // In addition, to maintain a manifold shell, we need to make sure that - // every edge from the input wire is used exactly once in the list of - // resulting faces. And that other internal edges are used twice. - - typedef std::pair uv_node; - - gp_Pln pln; - if (!approximate_plane_through_wire(wires.front(), pln, std::numeric_limits::infinity())) { - return false; - } - - const gp_XYZ& udir = pln.Position().XDirection().XYZ(); - const gp_XYZ& vdir = pln.Position().YDirection().XYZ(); - const gp_XYZ& pnt = pln.Position().Location().XYZ(); - - std::map mapping; - std::map, TopoDS_Edge> existing_edges, new_edges; - - std::unique_ptr mf; - - for (auto it = wires.begin(); it != wires.end(); ++it) { - const TopoDS_Wire& wire = *it; - BRepTools_WireExplorer exp(wire); - BRepBuilderAPI_MakePolygon mp; - - // Add UV coordinates to a newly created polygon - for (; exp.More(); exp.Next()) { - // Project onto plane - const TopoDS_Vertex& V = exp.CurrentVertex(); - gp_Pnt p = BRep_Tool::Pnt(V); - double u = (p.XYZ() - pnt).Dot(udir); - double v = (p.XYZ() - pnt).Dot(vdir); - mp.Add(gp_Pnt(u, v, 0.)); - - mapping.insert(std::make_pair(std::make_pair(u, v), V)); - - // Store existing edges in a map so that triangles can - // actually reference the preexisting edges. - const TopoDS_Edge& e = exp.Current(); - TopoDS_Vertex V0, V1; - TopExp::Vertices(e, V0, V1, true); - gp_Pnt p0 = BRep_Tool::Pnt(V0); - gp_Pnt p1 = BRep_Tool::Pnt(V1); - double u0 = (p0.XYZ() - pnt).Dot(udir); - double v0 = (p0.XYZ() - pnt).Dot(vdir); - double u1 = (p1.XYZ() - pnt).Dot(udir); - double v1 = (p1.XYZ() - pnt).Dot(vdir); - uv_node uv0 = std::make_pair(u0, v0); - uv_node uv1 = std::make_pair(u1, v1); - existing_edges.insert(std::make_pair(std::make_pair(uv0, uv1), e)); - existing_edges.insert(std::make_pair(std::make_pair(uv1, uv0), TopoDS::Edge(e.Reversed()))); - } - - // Not closed by default - mp.Close(); - - if (mf) { - if (it - 1 == wires.begin()) { - // @todo is this necessary? - TopoDS_Face f = mf->Face(); - mf->Init(f); - } - mf->Add(mp.Wire()); - } else { - mf.reset(new BRepBuilderAPI_MakeFace(mp.Wire())); - } - } - - const TopoDS_Face& face = mf->Face(); - - // Create a triangular mesh from the face - BRepMesh_IncrementalMesh(face, Precision::Confusion()); - - int n123[3]; - TopLoc_Location loc; - Handle_Poly_Triangulation tri = BRep_Tool::Triangulation(face, loc); - - if (!tri.IsNull()) { - - const Poly_Array1OfTriangle& triangles = tri->Triangles(); - for (int i = 1; i <= triangles.Length(); ++i) { - if (face.Orientation() == TopAbs_REVERSED) - triangles(i).Get(n123[2], n123[1], n123[0]); - else triangles(i).Get(n123[0], n123[1], n123[2]); - - // Create polygons from the mesh vertices - BRepBuilderAPI_MakeWire mp2; - for (int j = 0; j < 3; ++j) { - - uv_node uvnodes[2]; - TopoDS_Vertex vs[2]; - - for (int k = 0; k < 2; ++k) { - const gp_Pnt& uv = tri->Node(n123[(j + k) % 3]); - uvnodes[k] = std::make_pair(uv.X(), uv.Y()); - - auto it = mapping.find(uvnodes[k]); - if (it == mapping.end()) { - Logger::Error("Internal error: unable to unproject uv-mesh"); - return false; - } - - vs[k] = it->second; - } - - auto it = existing_edges.find(std::make_pair(uvnodes[0], uvnodes[1])); - if (it != existing_edges.end()) { - // This is a boundary edge, reuse existing edge from wire - mp2.Add(it->second); - } else { - auto jt = new_edges.find(std::make_pair(uvnodes[0], uvnodes[1])); - if (jt != new_edges.end()) { - // We have already added the reverse as part of another - // triangle, reuse this edge. - mp2.Add(TopoDS::Edge(jt->second)); - } else { - // This is a new internal edge. Register the reverse - // for reuse later. We need to be sure to reuse vertices - // for the edge construction because otherwise the wire - // builder will use geometrical proximity for vertex - // connections in which case the edge will be copied - // and no longer partner with other edges from the shell. - TopoDS_Edge ne = BRepBuilderAPI_MakeEdge(vs[0], vs[1]); - mp2.Add(ne); - // Store the reverse to be picked up later. - new_edges.insert(std::make_pair(std::make_pair(uvnodes[1], uvnodes[0]), TopoDS::Edge(ne.Reversed()))); - } - } - } - - BRepBuilderAPI_MakeFace mft(mp2.Wire()); - if (mft.IsDone()) { - TopoDS_Face triangle_face = mft.Face(); - TopoDS_Iterator jt(triangle_face, false); - for (; jt.More(); jt.Next()) { - const TopoDS_Wire& w = TopoDS::Wire(jt.Value()); - if (w.Orientation() != wires.front().Orientation()) { - triangle_face.Reverse(); - } - } - faces.Append(triangle_face); - } else { - Logger::Error("Internal error: missing face"); - return false; - } - } - } - - TopTools_IndexedDataMapOfShapeListOfShape mape, mapn; - for (auto& wire : wires) { - TopExp::MapShapesAndAncestors(wire, TopAbs_EDGE, TopAbs_WIRE, mape); - } - TopTools_ListIteratorOfListOfShape it(faces); - for (; it.More(); it.Next()) { - TopExp::MapShapesAndAncestors(it.Value(), TopAbs_EDGE, TopAbs_WIRE, mapn); - } - - // Validation - - for (int i = 1; i <= mape.Extent(); ++i) { -#if OCC_VERSION_HEX >= 0x70000 - TopTools_ListOfShape val; - if (!mapn.FindFromKey(mape.FindKey(i), val)) { -#else - bool contains = false; - try { - TopTools_ListOfShape val = mapn.FindFromKey(mape.FindKey(i)); - contains = true; - } catch (Standard_NoSuchObject&) {} - if (!contains) { -#endif - // All existing edges need to exist in the new faces - Logger::Error("Internal error, missing edge from triangulation"); - if (faceset_helper_ != nullptr) { - faceset_helper_->non_manifold() = true; - } - } - } - - for (int i = 1; i <= mapn.Extent(); ++i) { - const TopoDS_Shape& v = mapn.FindKey(i); - int n = mapn.FindFromIndex(i).Extent(); - // Existing edges are boundaries with use 1 - // New edges are internal with use 2 - if (n != (mape.Contains(v) ? 1 : 2)) { - Logger::Error("Internal error, non-manifold result from triangulation"); - if (faceset_helper_ != nullptr) { - faceset_helper_->non_manifold() = true; - } - } - } - - return true; -} - TopoDS_Shape IfcGeom::Kernel::apply_transformation(const TopoDS_Shape& s, const gp_Trsf& t) { if (t.Form() == gp_Identity) { return s; @@ -3597,260 +3312,6 @@ TopoDS_Shape IfcGeom::Kernel::apply_transformation(const TopoDS_Shape& s, const } } -namespace { - - /* - * A small helper utility to wrap around a numeric range - */ - class bounded_int { - private: - int i; - size_t n; - public: - bounded_int(int i, size_t n) : i(i), n(n) {} - - bounded_int& operator--() { - --i; - if (i == -1) { - i = (int) n - 1; - } - return *this; - } - - bounded_int& operator++() { - ++i; - if (i == (int) n) { - i = 0; - } - return *this; - } - - operator int() { return i; } - }; -} - -bool IfcGeom::Kernel::wire_intersections(const TopoDS_Wire& wire, TopTools_ListOfShape& wires) { - - if (getValue(GV_NO_WIRE_INTERSECTION_CHECK) > 0.) { - return false; - } - - if (!wire.Closed()) { - wires.Append(wire); - return false; - } - - int n = count(wire, TopAbs_EDGE); - if (n < 3) { - wires.Append(wire); - return false; - } - - // Note: initialize empty - Handle(ShapeExtend_WireData) wd = new ShapeExtend_WireData(); - - // ... to be sure to get consecutive edges - BRepTools_WireExplorer exp(wire); - IfcGeom::impl::tree tree; - - int edge_idx = 0; - for (; exp.More(); exp.Next()) { - wd->Add(exp.Current()); - if (n > 64) { - // tfk: indices in tree are 0-based vd 1-based in wiredata - tree.add(edge_idx++, exp.Current()); - } - } - - if (wd->NbEdges() != n) { - // If the number of edges differs, BRepTools_WireExplorer did not - // reach every edge, probably due to loops exactly at vertex locations. - // This is not supported by this algorithm which only elimates loops - // due to edge crossings. - - throw geometry_exception("Invalid loop"); - } - - bool intersected = false; - - // tfk: Extrema on infinite curves proved to be more robust. - // TopoDS_Face face = BRepBuilderAPI_MakeFace(wire, true).Face(); - // ShapeAnalysis_Wire saw(wd, face, getValue(GV_PRECISION)); - - double eps = 0; - if (getValue(GV_NO_WIRE_INTERSECTION_TOLERANCE) < 0.) { - eps = faceset_helper_ - // eps is added to both ends of the parametric domain, so 3. is chosen to be on the safe side here. - ? (faceset_helper_->epsilon() / 3.) - // @todo re-evaluate 2. here for the reasons above: - : (std::min)(util::min_edge_length(wire) / 2., getValue(GV_PRECISION) * 10.); - } - - // @todo: should this start from 0 in case of n > 64? - for (int i = 2; i < n; ++i) { - - std::vector js; - if (n > 64) { - Bnd_Box b; - BRepBndLib::Add(wd->Edge(i + 1), b); - b.Enlarge(eps); - js = tree.select_box(b, false); - } else { - boost::push_back(js, boost::irange(0, i - 1)); - } - - for(std::vector::const_iterator it = js.begin(); it != js.end(); ++it) { - int j = *it; - - if (n > 64) { - if (j > i) { - continue; - } - if ((std::max)(i, j) - (std::min)(i, j) <= 1) { - continue; - } - } - - // Only check non-consecutive edges - if (i == n - 1 && j == 0) continue; - - double u11, u12, u21, u22, U1, U2; - GeomAPI_ExtremaCurveCurve ecc( - BRep_Tool::Curve(wd->Edge(i + 1), u11, u12), - BRep_Tool::Curve(wd->Edge(j + 1), u21, u22) - ); - - // @todo: extend this to work in case of multiple extrema and curved segments. - const bool unbounded_intersects = (!ecc.Extrema().IsParallel() && ecc.NbExtrema() == 1 && ecc.Distance(1) < eps); - if (unbounded_intersects) { - ecc.Parameters(1, U1, U2); - - if (u11 > u12) { - std::swap(u11, u12); - } - if (u21 > u22) { - std::swap(u21, u22); - } - - /// @todo: tfk: probably need different thresholds on non-linear curves - u11 -= eps; - u12 += eps; - u21 -= eps; - u22 += eps; - - // tfk: code below is for ShapeAnalysis_Wire::CheckIntersectingEdges() - // IntRes2d_SequenceOfIntersectionPoint points2d; - // TColgp_SequenceOfPnt points3d; - // TColStd_SequenceOfReal errors; - // if (saw.CheckIntersectingEdges(i + 1, j + 1, points2d, points3d, errors)) { - - if (u11 < U1 && U1 < u12 && u21 < U2 && U2 < u22) { - - intersected = true; - - // Explore a forward and backward cycle from the intersection point - for (int fb = 0; fb <= 1; ++fb) { - const bool forward = fb == 0; - - BRepBuilderAPI_MakeWire mw; - bool first = true; - - for (bounded_int k(j, n);;) { - bool intersecting = k == j || k == i; - if (intersecting) { - TopoDS_Edge e = wd->Edge(k + 1); - - TopoDS_Vertex v1, v2; - TopExp::Vertices(e, v1, v2, true); - const TopoDS_Vertex* v = first == forward ? &v2 : &v1; - - // gp_Pnt p2 = points3d.Value(1); - - gp_Pnt p1 = BRep_Tool::Pnt(*v); - gp_Pnt pp1, pp2; - ecc.Points(1, pp1, pp2); - const gp_Pnt& p2 = k == i ? pp1 : pp2; - - // Substitute with a new edge from/to the intersection point - if (p1.Distance(p2) > getValue(GV_PRECISION) * 2) { - double _, __; - Handle_Geom_Curve crv = BRep_Tool::Curve(e, _, __); - BRepBuilderAPI_MakeEdge me(crv, p1, p2); - TopoDS_Edge ed = me.Edge(); - mw.Add(ed); - } - - first = false; - } else { - // Re-use original edge - mw.Add(wd->Edge(k + 1)); - } - - if (k == i) { - break; - } - - if (forward) { - ++k; - } else { - --k; - } - } - - // Recursively process both cuts - wire_intersections(mw.Wire(), wires); - } - - return true; - } - - } - } - } - - // No intersections found, append original wire - if (!intersected) { - wires.Append(wire); - } - - return intersected; -} - -void IfcGeom::Kernel::select_largest(const TopTools_ListOfShape& shapes, TopoDS_Shape& largest) { - double mass = 0.; - TopTools_ListIteratorOfListOfShape it(shapes); - for (; it.More(); it.Next()) { - /* - // tfk: bounding box is more efficient probably - const TopoDS_Wire& w = TopoDS::Wire(it.Value()); - TopoDS_Face face = BRepBuilderAPI_MakeFace(w).Face(); - const double m = face_area(face); - */ - - Bnd_Box bb; - BRepBndLib::AddClose(it.Value(), bb); - double xyz_min[3], xyz_max[3]; - bb.Get(xyz_min[0], xyz_min[1], xyz_min[2], xyz_max[0], xyz_max[1], xyz_max[2]); - const double eps = getValue(GV_PRECISION); - - double m = 1.; - for (int i = 0; i < 3; ++i) { - if (Precision::IsNegativeInfinite(xyz_min[i])) { - xyz_min[i] = 0.; - } - if (Precision::IsInfinite(xyz_max[i])) { - xyz_max[i] = 0.; - } - m *= (xyz_max[i] + eps) - (xyz_min[i] - eps); - } - - if (m > mass) { - mass = m; - largest = it.Value(); - } - } -} - bool IfcGeom::Kernel::fit_halfspace(const TopoDS_Shape& a, const TopoDS_Shape& b, TopoDS_Shape& box, double& height) { TopExp_Explorer exp(b, TopAbs_FACE); if (!exp.More()) { @@ -4105,6 +3566,16 @@ bool IfcGeom::Kernel::boolean_operation(const TopoDS_Shape& a_input, const TopTo return true; } + if (Logger::LOG_NOTICE >= Logger::Verbosity()) { + PERF("preliminary manifoldness check"); + + Logger::Notice("Operand A is " + (is_manifold(a) ? ""s : "non-"s) + "manifold"); + + TopTools_ListIteratorOfListOfShape it(b); + for (int i = 0; it.More(); it.Next(), ++i) { + Logger::Notice("Operand B " + std::to_string(i) + " is " + (is_manifold(it.Value()) ? ""s : "non-"s) + "manifold"); + } + } // Find a sensible value for the fuzziness, based on precision // and limited by edge lengths and vertex-edge distances. @@ -4543,361 +4014,6 @@ bool IfcGeom::Kernel::boolean_operation(const TopoDS_Shape& a, const TopoDS_Shap } #endif -namespace { - void find_neighbours(IfcGeom::impl::tree& tree, std::vector>& pnts, std::set& visited, int p, double eps) { - visited.insert(p); - - Bnd_Box b; - b.Set(*pnts[p].get()); - b.Enlarge(eps); - - std::vector js = tree.select_box(b, false); - for (int j : js) { - visited.insert(j); -#ifdef FACESET_HELPER_RECURSIVE - if (visited.find(j) == visited.end()) { - // @todo, making this recursive removes the dependence on the initial ordering, but will - // likely result in empty results when all vertices are within 1 eps from another point. - find_neighbours(tree, pnts, visited, j, eps); - } -#endif - } - } -} - -template -IfcGeom::Kernel::faceset_helper::~faceset_helper() { - // @todo this is super ugly, but how else can we be notified that the unique_ptr goes out of scope? - // Perhaps just supply a custom std::deleter? - kernel_->faceset_helper_ = nullptr; -} - - -template -bool IfcGeom::Kernel::faceset_helper::construct(const IfcSchema::IfcCartesianPoint* cp, gp_Pnt* l) { - return kernel_->convert(cp, *l); -} - -template -bool IfcGeom::Kernel::faceset_helper::construct(const std::vector& cp, gp_Pnt* l) { - if (cp.size() != 3) { - return false; - } - auto LU = kernel_->getValue(GV_LENGTH_UNIT); - l->SetCoord(cp[0] * LU, cp[1] * LU, cp[2] * LU); - return true; -} - -/* - -template -IfcGeom::Kernel::faceset_helper::faceset_helper(Kernel* kernel, const IfcSchema::IfcConnectedFaceSet* l) - : kernel_(kernel) - , non_manifold_(false) -{ - kernel->faceset_helper_ = this; - - IfcSchema::IfcCartesianPoint::list::ptr points = IfcParse::traverse((IfcUtil::IfcBaseClass*) l)->as(); - std::vector> pnts(std::distance(points->begin(), points->end())); - std::vector vertices(pnts.size()); - - IfcGeom::impl::tree tree; - - BRep_Builder B; - - Bnd_Box box; - for (size_t i = 0; i < points->size(); ++i) { - gp_Pnt* p = new gp_Pnt(); - if (kernel->convert(*(points->begin() + i), *p)) { - pnts[i].reset(p); - B.MakeVertex(vertices[i], *p, Precision::Confusion()); - tree.add(i, vertices[i]); - box.Add(*p); - } else { - delete p; - } - } - - // Use the bbox diagonal to influence local epsilon - // double bdiff = std::sqrt(box.SquareExtent()); - - // @todo the bounding box diagonal is not used (see above) - // because we're explicitly interested in the miminal - // dimension of the element to limit the tolerance (for sheet- - // like elements for example). But the way below is very - // dependent on orientation due to the usage of the - // axis-aligned bounding box. Use PCA to find three non-aligned - // set of dimensions and use the one with the smallest eigenvalue. - - // Find the minimal bounding box edge - double bmin[3], bmax[3]; - box.Get(bmin[0], bmin[1], bmin[2], bmax[0], bmax[1], bmax[2]); - double bdiff = std::numeric_limits::infinity(); - for (size_t i = 0; i < 3; ++i) { - const double d = bmax[i] - bmin[i]; - if (d > kernel->getValue(GV_PRECISION) * 10. && d < bdiff) { - bdiff = d; - } - } - - eps_ = kernel->getValue(GV_PRECISION) * 10. * (std::min)(1.0, bdiff); - - // @todo, there a tiny possibility that the duplicate faces are triggered - // for an internal boundary, that is also present as an external boundary. - // This will result in non-manifold configuration then, but this is deemed - // such as corner-case that it is not considered. - IfcSchema::IfcPolyLoop::list::ptr loops = IfcParse::traverse((IfcUtil::IfcBaseClass*)l)->as(); - - size_t loops_removed, non_manifold, duplicate_faces; - - std::map, int> edge_use; - - for (int i = 0; i < 3; ++i) { - // Some times files, have large tolerance values specified collapsing too many vertices. - // This case we detect below and re-run the loop with smaller epsilon. Normally - // the body of this loop would only be executed once. - - loops_removed = 0; - non_manifold = 0; - duplicate_faces = 0; - - vertex_mapping_.clear(); - duplicates_.clear(); - - edge_use.clear(); - - if (eps_ < Precision::Confusion()) { - // occt uses some hard coded precision values, don't go smaller than that. - // @todo, can be reset though with BRepLib::Precision(double) - eps_ = Precision::Confusion(); - } - - for (int pnt_i = 0; pnt_i < (int)pnts.size(); ++pnt_i) { - if (pnts[pnt_i]) { - std::set vs; - find_neighbours(tree, pnts, vs, pnt_i, eps_); - - for (int v : vs) { - auto pt = *(points->begin() + v); - // NB: insert() ignores duplicate keys - vertex_mapping_.insert({ get_idx(pt), pnt_i }); - } - } - } - - typedef std::array edge_t; - typedef std::set edge_set_t; - std::set edge_sets; - - for (auto& loop : *loops) { - auto ps = loop->Polygon(); - - std::vector > segments; - edge_set_t segment_set; - - loop_(ps, [&segments, &segment_set](int C, int D, bool) { - segment_set.insert(edge_t{C,D}); - segments.push_back(std::make_pair(C, D)); - }); - - if (edge_sets.find(segment_set) != edge_sets.end()) { - duplicate_faces++; - duplicates_.insert(loop); - continue; - } - edge_sets.insert(segment_set); - - if (segments.size() >= 3) { - for (auto& p : segments) { - edge_use[p] ++; - } - } else { - loops_removed += 1; - } - } - - if (edge_use.size() != 0) { - break; - } else { - eps_ /= 10.; - } - } - - for (auto& p : edge_use) { - int a, b; - std::tie(a, b) = p.first; - edges_[p.first] = BRepBuilderAPI_MakeEdge(vertices[a], vertices[b]); - - if (p.second != 2) { - non_manifold += 1; - } - } - - if (loops_removed || (non_manifold && l->declaration().is(IfcSchema::IfcClosedShell::Class()))) { - Logger::Warning(boost::lexical_cast(duplicate_faces) + " duplicate faces removed, " + boost::lexical_cast(loops_removed) + " loops removed and " + boost::lexical_cast(non_manifold) + " non-manifold edges for:", l); - } -} - -*/ - -namespace { - const std::vector>* store_cache(const std::vector>& p) { - return &p; - } - - const std::vector>* store_cache(const std::vector& /*p*/) { - return nullptr; - } -} - -template -IfcGeom::Kernel::faceset_helper::faceset_helper( - Kernel* kernel, - const std::vector& points, - const std::vector& indices, - bool should_be_closed -) - : kernel_(kernel) - , non_manifold_(false) - , points_(store_cache(points)) -{ - std::vector> pnts(std::distance(points.begin(), points.end())); - std::vector vertices(pnts.size()); - - IfcGeom::impl::tree tree; - - BRep_Builder B; - - Bnd_Box box; - for (size_t i = 0; i < points.size(); ++i) { - gp_Pnt* p = new gp_Pnt; - if (construct(points[i], p)) { - pnts[i].reset(p); - B.MakeVertex(vertices[i], *p, Precision::Confusion()); - tree.add((int) i, vertices[i]); - box.Add(*p); - } else { - delete p; - } - } - - // Use the bbox diagonal to influence local epsilon - // double bdiff = std::sqrt(box.SquareExtent()); - - // @todo the bounding box diagonal is not used (see above) - // because we're explicitly interested in the miminal - // dimension of the element to limit the tolerance (for sheet- - // like elements for example). But the way below is very - // dependent on orientation due to the usage of the - // axis-aligned bounding box. Use PCA to find three non-aligned - // set of dimensions and use the one with the smallest eigenvalue. - - // Find the minimal bounding box edge - double bmin[3], bmax[3]; - box.Get(bmin[0], bmin[1], bmin[2], bmax[0], bmax[1], bmax[2]); - double bdiff = std::numeric_limits::infinity(); - for (size_t i = 0; i < 3; ++i) { - const double d = bmax[i] - bmin[i]; - if (d > kernel->getValue(GV_PRECISION) * 10. && d < bdiff) { - bdiff = d; - } - } - - eps_ = kernel->getValue(GV_PRECISION) * 10. * (std::min)(1.0, bdiff); - - size_t loops_removed, non_manifold, duplicate_faces; - - std::map, int> edge_use; - - for (int i = 0; i < 3; ++i) { - // Some times files, have large tolerance values specified collapsing too many vertices. - // This case we detect below and re-run the loop with smaller epsilon. Normally - // the body of this loop would only be executed once. - - loops_removed = 0; - non_manifold = 0; - duplicate_faces = 0; - - vertex_mapping_.clear(); - duplicates_.clear(); - - edge_use.clear(); - - if (eps_ < Precision::Confusion()) { - // occt uses some hard coded precision values, don't go smaller than that. - // @todo, can be reset though with BRepLib::Precision(double) - eps_ = Precision::Confusion(); - } - - for (int pnt_i = 0; pnt_i < (int)pnts.size(); ++pnt_i) { - if (pnts[pnt_i]) { - std::set vs; - find_neighbours(tree, pnts, vs, pnt_i, eps_); - - for (int v : vs) { - // NB: insert() ignores duplicate keys - // v-1? - vertex_mapping_.insert({ get_idx(points[v]), pnt_i }); - } - } - } - - typedef std::array edge_t; - typedef std::set edge_set_t; - std::set edge_sets; - - for (auto ps = indices.begin(); ps != indices.end(); ++ps) { - std::vector > segments; - edge_set_t segment_set; - - loop_(*ps, [&segments, &segment_set](int C, int D, bool) { - segment_set.insert(edge_t{ C,D }); - segments.push_back(std::make_pair(C, D)); - }); - - if (edge_sets.find(segment_set) != edge_sets.end()) { - duplicate_faces++; - duplicates_.insert(util::conditional_address_of(*ps)); - continue; - } - edge_sets.insert(segment_set); - - if (segments.size() >= 3) { - for (auto& p : segments) { - edge_use[p] ++; - } - } - else { - loops_removed += 1; - } - } - - if (edge_use.size() != 0) { - break; - } - else { - eps_ /= 10.; - } - } - - for (auto& p : edge_use) { - int a, b; - std::tie(a, b) = p.first; - edges_[p.first] = BRepBuilderAPI_MakeEdge(vertices[a], vertices[b]); - - if (p.second != 2) { - non_manifold += 1; - } - } - - if (duplicates_.size() || loops_removed || (non_manifold && should_be_closed)) { - Logger::Warning(boost::lexical_cast(duplicate_faces) + " duplicate faces removed, " + boost::lexical_cast(loops_removed) + " degenerate loops eliminated and " + boost::lexical_cast(non_manifold) + " non-manifold edges"); - } -} - -template class IfcGeom::Kernel::faceset_helper; -template class IfcGeom::Kernel::faceset_helper, std::vector>; - void IfcGeom::Kernel::set_conversion_placement_rel_to_type(const IfcParse::declaration* type) { placement_rel_to_type_ = type; } diff --git a/src/ifcgeom/IfcGeom.h b/src/ifcgeom/IfcGeom.h index ae0ec58d06..d9a79353ec 100644 --- a/src/ifcgeom/IfcGeom.h +++ b/src/ifcgeom/IfcGeom.h @@ -119,28 +119,8 @@ private: const std::vector>* points_ = nullptr; double eps_; bool non_manifold_; - - template - void loop_(const LP& lp, const Fn& callback) { - auto ps = get_idxs(lp); - - if (ps.size() < 3) { - return; - } - - auto A = ps.back(); - for (auto& B : ps) { - auto C = vertex_mapping_[A], D = vertex_mapping_[B]; - bool fwd = C < D; - if (!fwd) { - std::swap(C, D); - } - if (C != D) { - callback(C, D, fwd); - A = B; - } - } - } + + void loop_(const LP& lp, const std::function& callback); bool construct(const IfcSchema::IfcCartesianPoint* cp, gp_Pnt* l); bool construct(const std::vector& cp, gp_Pnt* l); @@ -153,32 +133,9 @@ private: return &cp; } - std::vector get_idxs(const IfcSchema::IfcPolyLoop* lp) { - auto poly = lp->Polygon(); - std::vector idxs; - std::transform(poly->begin(), poly->end(), std::back_inserter(idxs), [this](const IfcSchema::IfcCartesianPoint* p) {return get_idx(p); }); - return idxs; - } - - std::vector get_idxs(const std::vector& it) { - std::vector idxs; - std::transform(it.begin(), it.end(), std::back_inserter(idxs), [this](int i) { return get_idx((*points_)[i - 1]); }); - return idxs; - } - - /* - std::vector get_idxs(std::vector>::const_iterator it) { - std::vector idxs; - std::transform(it->begin(), it->end(), std::back_inserter(idxs), [this](int i) {return get_idx(i); }); - return idxs; - } - */ - + std::vector get_idxs(const IfcSchema::IfcPolyLoop* lp); + std::vector get_idxs(const std::vector& it); public: - /* - faceset_helper(MAKE_TYPE_NAME(Kernel)* kernel, const IfcSchema::IfcConnectedFaceSet* l); - */ - faceset_helper( MAKE_TYPE_NAME(Kernel)* kernel, const std::vector& points, @@ -189,64 +146,11 @@ private: bool non_manifold() const { return non_manifold_; } bool& non_manifold() { return non_manifold_; } + double epsilon() const { return eps_; } + + bool edge(int A, int B, TopoDS_Edge& e); - /* - bool edge(CP a, CP b, TopoDS_Edge& e) { - int A = vertex_mapping_[get_idx(a)]; - int B = vertex_mapping_[get_idx(b)]; - if (A == B) { - return false; - } - - return edge(A, B, e); - } - */ - - bool edge(int A, int B, TopoDS_Edge& e) { - auto it = edges_.find({A, B}); - if (it == edges_.end()) { - return false; - } - e = it->second; - return true; - } - - bool wire(const LP& loop, TopoDS_Wire& wire) { - if (duplicates_.find(util::conditional_address_of(loop)) != duplicates_.end()) { - return false; - } - BRep_Builder builder; - builder.MakeWire(wire); - int count = 0; - loop_(loop, [this, &builder, &wire, &count](int A, int B, bool fwd) { - TopoDS_Edge e; - if (edge(A, B, e)) { - if (!fwd) { - e.Reverse(); - } - builder.Add(wire, e); - count += 1; - } - }); - if (count >= 3) { - wire.Closed(true); - - TopTools_ListOfShape results; - if (kernel_->wire_intersections(wire, results)) { - Logger::Warning("Self-intersections with " + boost::lexical_cast(results.Extent()) + " cycles detected"); - kernel_->select_largest(results, wire); - non_manifold_ = true; - } - - return true; - } else { - return false; - } - } - - double epsilon() const { - return eps_; - } + bool wire(const LP& loop, TopoDS_Wire& wire); }; double deflection_tolerance; @@ -355,6 +259,7 @@ public: void set_offset(const std::array& offset); void set_rotation(const std::array& rotation); + double get_wire_intersection_tolerance(const TopoDS_Wire&) const; bool convert_wire_to_face(const TopoDS_Wire& wire, TopoDS_Face& face); bool convert_wire_to_faces(const TopoDS_Wire& wire, TopoDS_Compound& face); @@ -419,12 +324,6 @@ public: void remove_collinear_points_from_loop(TColgp_SequenceOfPnt& polygon, bool closed, double tol=-1.); bool wire_to_sequence_of_point(const TopoDS_Wire&, TColgp_SequenceOfPnt&); void sequence_of_point_to_wire(const TColgp_SequenceOfPnt&, TopoDS_Wire&, bool closed); - bool approximate_plane_through_wire(const TopoDS_Wire&, gp_Pln&, double eps=-1.); - bool flatten_wire(TopoDS_Wire&); - /// Triangulate the set of wires. The firstmost wire is assumed to be the outer wire. - bool triangulate_wire(const std::vector&, TopTools_ListOfShape&); - bool wire_intersections(const TopoDS_Wire & wire, TopTools_ListOfShape & wires); - void select_largest(const TopTools_ListOfShape& shapes, TopoDS_Shape& largest); static double shape_volume(const TopoDS_Shape& s); static double face_area(const TopoDS_Face& f); diff --git a/src/ifcgeom/IfcPolyLoop.cpp b/src/ifcgeom/IfcPolyLoop.cpp index 7a1b5b4deb..49a0b3a683 100644 --- a/src/ifcgeom/IfcPolyLoop.cpp +++ b/src/ifcgeom/IfcPolyLoop.cpp @@ -17,11 +17,13 @@ * * ********************************************************************************/ +#include "../ifcgeom/IfcGeom.h" +#include "../ifcgeom_schema_agnostic/wire_utils.h" + #include #include #include #include -#include "../ifcgeom/IfcGeom.h" #define _USE_MATH_DEFINES #define Kernel MAKE_TYPE_NAME(Kernel) @@ -68,9 +70,9 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcPolyLoop* l, TopoDS_Wire& resu result = w.Wire(); TopTools_ListOfShape results; - if (wire_intersections(result, results)) { + if (getValue(GV_NO_WIRE_INTERSECTION_CHECK) == 0. && util::wire_intersections(result, results, get_wire_intersection_tolerance(result), getValue(GV_PRECISION))) { Logger::Error("Self-intersections with " + boost::lexical_cast(results.Extent()) + " cycles detected", l); - select_largest(results, result); + util::select_largest(results, result); } return true; diff --git a/src/ifcgeom/IfcPolygonalFaceSet.cpp b/src/ifcgeom/IfcPolygonalFaceSet.cpp index 1ef9e99c2b..4703a11df5 100644 --- a/src/ifcgeom/IfcPolygonalFaceSet.cpp +++ b/src/ifcgeom/IfcPolygonalFaceSet.cpp @@ -22,6 +22,7 @@ #include #include #include "../ifcgeom/IfcGeom.h" +#include "../ifcgeom_schema_agnostic/wire_utils.h" #define Kernel MAKE_TYPE_NAME(Kernel) @@ -111,7 +112,7 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcPolygonalFaceSet* pfs, TopoDS_ if (not_planar) { TopTools_ListOfShape fs; - if (triangulate_wire(ws, fs)) { + if (util::triangulate_wire(ws, fs)) { Logger::Warning("Triangulated face boundary:", pfs); TopTools_ListIteratorOfListOfShape it(fs); for (; it.More(); it.Next()) { diff --git a/src/ifcgeom/faceset_helper.cpp b/src/ifcgeom/faceset_helper.cpp new file mode 100644 index 0000000000..62c6ea747c --- /dev/null +++ b/src/ifcgeom/faceset_helper.cpp @@ -0,0 +1,301 @@ +#include "IfcGeom.h" + +#include "../ifcgeom_schema_agnostic/IfcGeomTree.h" +#include "../ifcgeom_schema_agnostic/wire_utils.h" + +#define Kernel MAKE_TYPE_NAME(Kernel) + +namespace { + void find_neighbours(IfcGeom::impl::tree& tree, std::vector>& pnts, std::set& visited, int p, double eps) { + visited.insert(p); + + Bnd_Box b; + b.Set(*pnts[p].get()); + b.Enlarge(eps); + + std::vector js = tree.select_box(b, false); + for (int j : js) { + visited.insert(j); +#ifdef FACESET_HELPER_RECURSIVE + if (visited.find(j) == visited.end()) { + // @todo, making this recursive removes the dependence on the initial ordering, but will + // likely result in empty results when all vertices are within 1 eps from another point. + find_neighbours(tree, pnts, visited, j, eps); + } +#endif + } + } +} + +namespace { + const std::vector>* store_cache(const std::vector>& p) { + return &p; + } + + const std::vector>* store_cache(const std::vector& /*p*/) { + return nullptr; + } +} + +template +IfcGeom::Kernel::faceset_helper::faceset_helper( + Kernel* kernel, + const std::vector& points, + const std::vector& indices, + bool should_be_closed +) + : kernel_(kernel) + , non_manifold_(false) + , points_(store_cache(points)) +{ + std::vector> pnts(std::distance(points.begin(), points.end())); + std::vector vertices(pnts.size()); + + IfcGeom::impl::tree tree; + + BRep_Builder B; + + Bnd_Box box; + for (size_t i = 0; i < points.size(); ++i) { + gp_Pnt* p = new gp_Pnt; + if (construct(points[i], p)) { + pnts[i].reset(p); + B.MakeVertex(vertices[i], *p, Precision::Confusion()); + tree.add((int)i, vertices[i]); + box.Add(*p); + } else { + delete p; + } + } + + // Use the bbox diagonal to influence local epsilon + // double bdiff = std::sqrt(box.SquareExtent()); + + // @todo the bounding box diagonal is not used (see above) + // because we're explicitly interested in the miminal + // dimension of the element to limit the tolerance (for sheet- + // like elements for example). But the way below is very + // dependent on orientation due to the usage of the + // axis-aligned bounding box. Use PCA to find three non-aligned + // set of dimensions and use the one with the smallest eigenvalue. + + // Find the minimal bounding box edge + double bmin[3], bmax[3]; + box.Get(bmin[0], bmin[1], bmin[2], bmax[0], bmax[1], bmax[2]); + double bdiff = std::numeric_limits::infinity(); + for (size_t i = 0; i < 3; ++i) { + const double d = bmax[i] - bmin[i]; + if (d > kernel->getValue(GV_PRECISION) * 10. && d < bdiff) { + bdiff = d; + } + } + + eps_ = kernel->getValue(GV_PRECISION) * 10. * (std::min)(1.0, bdiff); + + size_t loops_removed, non_manifold, duplicate_faces; + + std::map, int> edge_use; + + for (int i = 0; i < 3; ++i) { + // Some times files, have large tolerance values specified collapsing too many vertices. + // This case we detect below and re-run the loop with smaller epsilon. Normally + // the body of this loop would only be executed once. + + loops_removed = 0; + non_manifold = 0; + duplicate_faces = 0; + + vertex_mapping_.clear(); + duplicates_.clear(); + + edge_use.clear(); + + if (eps_ < Precision::Confusion()) { + // occt uses some hard coded precision values, don't go smaller than that. + // @todo, can be reset though with BRepLib::Precision(double) + eps_ = Precision::Confusion(); + } + + for (int pnt_i = 0; pnt_i < (int)pnts.size(); ++pnt_i) { + if (pnts[pnt_i]) { + std::set vs; + find_neighbours(tree, pnts, vs, pnt_i, eps_); + + for (int v : vs) { + // NB: insert() ignores duplicate keys + // v-1? + vertex_mapping_.insert({ get_idx(points[v]), pnt_i }); + } + } + } + + std::set> unique; + for (int pnt_i = 0; pnt_i < (int)pnts.size(); ++pnt_i) { + if (pnts[pnt_i]) { + unique.insert(std::make_tuple( + (*pnts[pnt_i]).X(), + (*pnts[pnt_i]).Y(), + (*pnts[pnt_i]).Z() + )); + } + } + + Logger::Notice("Collapsed vertices from " + std::to_string(pnts.size()) + " (" + std::to_string(unique.size()) + " unique) to " + std::to_string(vertex_mapping_.size())); + + typedef std::array edge_t; + typedef std::set edge_set_t; + std::set edge_sets; + + for (auto ps = indices.begin(); ps != indices.end(); ++ps) { + std::vector > segments; + edge_set_t segment_set; + + loop_(*ps, [&segments, &segment_set](int C, int D, bool) { + segment_set.insert(edge_t{ C,D }); + segments.push_back(std::make_pair(C, D)); + }); + + if (edge_sets.find(segment_set) != edge_sets.end()) { + duplicate_faces++; + duplicates_.insert(util::conditional_address_of(*ps)); + continue; + } + edge_sets.insert(segment_set); + + if (segments.size() >= 3) { + for (auto& p : segments) { + edge_use[p] ++; + } + } else { + loops_removed += 1; + } + } + + if (edge_use.size() != 0) { + break; + } else { + eps_ /= 10.; + } + } + + for (auto& p : edge_use) { + int a, b; + std::tie(a, b) = p.first; + edges_[p.first] = BRepBuilderAPI_MakeEdge(vertices[a], vertices[b]); + + if (p.second != 2) { + non_manifold += 1; + } + } + + if (duplicates_.size() || loops_removed || (non_manifold && should_be_closed)) { + Logger::Warning(boost::lexical_cast(duplicate_faces) + " duplicate faces removed, " + boost::lexical_cast(loops_removed) + " degenerate loops eliminated and " + boost::lexical_cast(non_manifold) + " non-manifold edges"); + } +} + +template +void IfcGeom::Kernel::faceset_helper::loop_(const LP& lp, const std::function& callback) { + auto ps = get_idxs(lp); + + if (ps.size() < 3) { + return; + } + + auto A = ps.back(); + for (auto& B : ps) { + auto C = vertex_mapping_[A], D = vertex_mapping_[B]; + bool fwd = C < D; + if (!fwd) { + std::swap(C, D); + } + if (C != D) { + callback(C, D, fwd); + A = B; + } + } +} + +template +std::vector IfcGeom::Kernel::faceset_helper::get_idxs(const IfcSchema::IfcPolyLoop* lp) { + auto poly = lp->Polygon(); + std::vector idxs; + std::transform(poly->begin(), poly->end(), std::back_inserter(idxs), [this](const IfcSchema::IfcCartesianPoint* p) {return get_idx(p); }); + return idxs; +} + +template +std::vector IfcGeom::Kernel::faceset_helper::get_idxs(const std::vector& it) { + std::vector idxs; + std::transform(it.begin(), it.end(), std::back_inserter(idxs), [this](int i) { return get_idx((*points_)[i - 1]); }); + return idxs; +} + +template +bool IfcGeom::Kernel::faceset_helper::edge(int A, int B, TopoDS_Edge& e) { + auto it = edges_.find({ A, B }); + if (it == edges_.end()) { + return false; + } + e = it->second; + return true; +} + +template +bool IfcGeom::Kernel::faceset_helper::wire(const LP& loop, TopoDS_Wire& wire) { + if (duplicates_.find(util::conditional_address_of(loop)) != duplicates_.end()) { + return false; + } + BRep_Builder builder; + builder.MakeWire(wire); + int count = 0; + loop_(loop, [this, &builder, &wire, &count](int A, int B, bool fwd) { + TopoDS_Edge e; + if (edge(A, B, e)) { + if (!fwd) { + e.Reverse(); + } + builder.Add(wire, e); + count += 1; + } + }); + if (count >= 3) { + wire.Closed(true); + + TopTools_ListOfShape results; + if (kernel_->getValue(GV_NO_WIRE_INTERSECTION_CHECK) == 0. && util::wire_intersections(wire, results, kernel_->get_wire_intersection_tolerance(wire), kernel_->getValue(IfcGeom::Kernel::GV_PRECISION))) { + Logger::Warning("Self-intersections with " + boost::lexical_cast(results.Extent()) + " cycles detected"); + util::select_largest(results, wire); + non_manifold_ = true; + } + + return true; + } else { + return false; + } +} + +template +IfcGeom::Kernel::faceset_helper::~faceset_helper() { + // @todo this is super ugly, but how else can we be notified that the unique_ptr goes out of scope? + // Perhaps just supply a custom std::deleter? + kernel_->faceset_helper_ = nullptr; +} + + +template +bool IfcGeom::Kernel::faceset_helper::construct(const IfcSchema::IfcCartesianPoint* cp, gp_Pnt* l) { + return kernel_->convert(cp, *l); +} + +template +bool IfcGeom::Kernel::faceset_helper::construct(const std::vector& cp, gp_Pnt* l) { + if (cp.size() != 3) { + return false; + } + auto LU = kernel_->getValue(GV_LENGTH_UNIT); + l->SetCoord(cp[0] * LU, cp[1] * LU, cp[2] * LU); + return true; +} + +template class IfcGeom::Kernel::faceset_helper; +template class IfcGeom::Kernel::faceset_helper, std::vector>; \ No newline at end of file diff --git a/src/ifcgeom_schema_agnostic/wire_utils.cpp b/src/ifcgeom_schema_agnostic/wire_utils.cpp new file mode 100644 index 0000000000..2c069bbe27 --- /dev/null +++ b/src/ifcgeom_schema_agnostic/wire_utils.cpp @@ -0,0 +1,558 @@ +#include "wire_utils.h" + +#include "../ifcparse/IfcLogger.h" +#include "../ifcgeom_schema_agnostic/Kernel.h" +#include "../ifcgeom_schema_agnostic/IfcGeomTree.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +bool IfcGeom::util::approximate_plane_through_wire(const TopoDS_Wire& wire, gp_Pln& plane, double eps_) { + // Newell's Method is used for the normal calculation + // as a simple edge cross product can give opposite results + // for a concave face boundary. + // Reference: Graphics Gems III p. 231 + + const double eps2 = eps_ * eps_; + + double x = 0, y = 0, z = 0; + gp_Pnt current, previous, first; + gp_XYZ center; + int n = 0; + + BRepTools_WireExplorer exp(wire); + + for (;; exp.Next()) { + const bool has_more = exp.More() != 0; + if (has_more) { + const TopoDS_Vertex& v = exp.CurrentVertex(); + current = BRep_Tool::Pnt(v); + center += current.XYZ(); + } else { + current = first; + } + if (n) { + const double& xn = previous.X(); + const double& yn = previous.Y(); + const double& zn = previous.Z(); + const double& xn1 = current.X(); + const double& yn1 = current.Y(); + const double& zn1 = current.Z(); + x += (yn - yn1)*(zn + zn1); + y += (xn + xn1)*(zn - zn1); + z += (xn - xn1)*(yn + yn1); + } else { + first = current; + } + if (!has_more) { + break; + } + previous = current; + ++n; + } + + if (n < 3) { + return false; + } + + plane = gp_Pln(center / n, gp_Dir(x, y, z)); + + exp.Init(wire); + for (; exp.More(); exp.Next()) { + const TopoDS_Vertex& v = exp.CurrentVertex(); + current = BRep_Tool::Pnt(v); + if (plane.SquareDistance(current) > eps2) { + return false; + } + } + + return true; +} + +bool IfcGeom::util::flatten_wire(TopoDS_Wire& wire, double eps) { + gp_Pln pln; + if (!approximate_plane_through_wire(wire, pln, eps)) { + return false; + } + TopoDS_Face face = BRepBuilderAPI_MakeFace(pln).Face(); + BRepAlgo_NormalProjection proj(face); + proj.Add(wire); + proj.Build(); + if (!proj.IsDone()) { + return false; + } + TopTools_ListOfShape list; + proj.BuildWire(list); + if (list.Extent() != 1) { + return false; + } + wire = TopoDS::Wire(list.First()); + return true; +} + +IfcGeom::util::triangulate_wire_result IfcGeom::util::triangulate_wire(const std::vector& wires, TopTools_ListOfShape& faces) { + // This is a bit of a precarious approach, but seems to work for the + // versions of OCCT tested for. OCCT has a Delaunay triangulation function + // BRepMesh_Delaun, but it is notoriously hard to interpret the results + // (due to the Bowyer-Watson super triangle perhaps?). Therefore + // alternatively we use the regular OCCT incremental mesher on a new face + // created from the UV coordinates of the original wire. Pray to our gods + // that the vertex coordinates are unaffected by the meshing algorithm and + // map them back to 3d coordinates when iterating over the mesh triangles. + + // In addition, to maintain a manifold shell, we need to make sure that + // every edge from the input wire is used exactly once in the list of + // resulting faces. And that other internal edges are used twice. + + typedef std::pair uv_node; + + gp_Pln pln; + if (!approximate_plane_through_wire(wires.front(), pln, std::numeric_limits::infinity())) { + return TRIANGULATE_WIRE_FAIL; + } + + const gp_XYZ& udir = pln.Position().XDirection().XYZ(); + const gp_XYZ& vdir = pln.Position().YDirection().XYZ(); + const gp_XYZ& pnt = pln.Position().Location().XYZ(); + + std::map mapping; + std::map, TopoDS_Edge> existing_edges, new_edges; + + std::unique_ptr mf; + + for (auto it = wires.begin(); it != wires.end(); ++it) { + const TopoDS_Wire& wire = *it; + BRepTools_WireExplorer exp(wire); + BRepBuilderAPI_MakePolygon mp; + + // Add UV coordinates to a newly created polygon + for (; exp.More(); exp.Next()) { + // Project onto plane + const TopoDS_Vertex& V = exp.CurrentVertex(); + gp_Pnt p = BRep_Tool::Pnt(V); + double u = (p.XYZ() - pnt).Dot(udir); + double v = (p.XYZ() - pnt).Dot(vdir); + mp.Add(gp_Pnt(u, v, 0.)); + + mapping.insert(std::make_pair(std::make_pair(u, v), V)); + + // Store existing edges in a map so that triangles can + // actually reference the preexisting edges. + const TopoDS_Edge& e = exp.Current(); + TopoDS_Vertex V0, V1; + TopExp::Vertices(e, V0, V1, true); + gp_Pnt p0 = BRep_Tool::Pnt(V0); + gp_Pnt p1 = BRep_Tool::Pnt(V1); + double u0 = (p0.XYZ() - pnt).Dot(udir); + double v0 = (p0.XYZ() - pnt).Dot(vdir); + double u1 = (p1.XYZ() - pnt).Dot(udir); + double v1 = (p1.XYZ() - pnt).Dot(vdir); + uv_node uv0 = std::make_pair(u0, v0); + uv_node uv1 = std::make_pair(u1, v1); + existing_edges.insert(std::make_pair(std::make_pair(uv0, uv1), e)); + existing_edges.insert(std::make_pair(std::make_pair(uv1, uv0), TopoDS::Edge(e.Reversed()))); + } + + // Not closed by default + mp.Close(); + + if (mf) { + if (it - 1 == wires.begin()) { + // @todo is this necessary? + TopoDS_Face f = mf->Face(); + mf->Init(f); + } + mf->Add(mp.Wire()); + } else { + mf.reset(new BRepBuilderAPI_MakeFace(mp.Wire())); + } + } + + const TopoDS_Face& face = mf->Face(); + + // Create a triangular mesh from the face + BRepMesh_IncrementalMesh(face, Precision::Confusion()); + + int n123[3]; + TopLoc_Location loc; + Handle_Poly_Triangulation tri = BRep_Tool::Triangulation(face, loc); + + if (!tri.IsNull()) { + + const Poly_Array1OfTriangle& triangles = tri->Triangles(); + for (int i = 1; i <= triangles.Length(); ++i) { + if (face.Orientation() == TopAbs_REVERSED) + triangles(i).Get(n123[2], n123[1], n123[0]); + else triangles(i).Get(n123[0], n123[1], n123[2]); + + // Create polygons from the mesh vertices + BRepBuilderAPI_MakeWire mp2; + for (int j = 0; j < 3; ++j) { + + uv_node uvnodes[2]; + TopoDS_Vertex vs[2]; + + for (int k = 0; k < 2; ++k) { + const gp_Pnt& uv = tri->Node(n123[(j + k) % 3]); + uvnodes[k] = std::make_pair(uv.X(), uv.Y()); + + auto it = mapping.find(uvnodes[k]); + if (it == mapping.end()) { + Logger::Error("Internal error: unable to unproject uv-mesh"); + return TRIANGULATE_WIRE_FAIL; + } + + vs[k] = it->second; + } + + auto it = existing_edges.find(std::make_pair(uvnodes[0], uvnodes[1])); + if (it != existing_edges.end()) { + // This is a boundary edge, reuse existing edge from wire + mp2.Add(it->second); + } else { + auto jt = new_edges.find(std::make_pair(uvnodes[0], uvnodes[1])); + if (jt != new_edges.end()) { + // We have already added the reverse as part of another + // triangle, reuse this edge. + mp2.Add(TopoDS::Edge(jt->second)); + } else { + // This is a new internal edge. Register the reverse + // for reuse later. We need to be sure to reuse vertices + // for the edge construction because otherwise the wire + // builder will use geometrical proximity for vertex + // connections in which case the edge will be copied + // and no longer partner with other edges from the shell. + TopoDS_Edge ne = BRepBuilderAPI_MakeEdge(vs[0], vs[1]); + mp2.Add(ne); + // Store the reverse to be picked up later. + new_edges.insert(std::make_pair(std::make_pair(uvnodes[1], uvnodes[0]), TopoDS::Edge(ne.Reversed()))); + } + } + } + + BRepBuilderAPI_MakeFace mft(mp2.Wire()); + if (mft.IsDone()) { + TopoDS_Face triangle_face = mft.Face(); + TopoDS_Iterator jt(triangle_face, false); + for (; jt.More(); jt.Next()) { + const TopoDS_Wire& w = TopoDS::Wire(jt.Value()); + if (w.Orientation() != wires.front().Orientation()) { + triangle_face.Reverse(); + } + } + faces.Append(triangle_face); + } else { + Logger::Error("Internal error: missing face"); + return TRIANGULATE_WIRE_FAIL; + } + } + } + + TopTools_IndexedDataMapOfShapeListOfShape mape, mapn; + for (auto& wire : wires) { + TopExp::MapShapesAndAncestors(wire, TopAbs_EDGE, TopAbs_WIRE, mape); + } + TopTools_ListIteratorOfListOfShape it(faces); + for (; it.More(); it.Next()) { + TopExp::MapShapesAndAncestors(it.Value(), TopAbs_EDGE, TopAbs_WIRE, mapn); + } + + // Validation + bool non_manifold = false; + + for (int i = 1; i <= mape.Extent(); ++i) { +#if OCC_VERSION_HEX >= 0x70000 + TopTools_ListOfShape val; + if (!mapn.FindFromKey(mape.FindKey(i), val)) { +#else + bool contains = false; + try { + TopTools_ListOfShape val = mapn.FindFromKey(mape.FindKey(i)); + contains = true; + } catch (Standard_NoSuchObject&) {} + if (!contains) { +#endif + // All existing edges need to exist in the new faces + Logger::Error("Internal error, missing edge from triangulation"); + non_manifold = true; + } + } + + for (int i = 1; i <= mapn.Extent(); ++i) { + const TopoDS_Shape& v = mapn.FindKey(i); + int n = mapn.FindFromIndex(i).Extent(); + // Existing edges are boundaries with use 1 + // New edges are internal with use 2 + if (n != (mape.Contains(v) ? 1 : 2)) { + Logger::Error("Internal error, non-manifold result from triangulation"); + non_manifold = true; + } + } + + return non_manifold ? TRIANGULATE_WIRE_NON_MANIFOLD : TRIANGULATE_WIRE_OK; +} + +namespace { + + /* + * A small helper utility to wrap around a numeric range + */ + class bounded_int { + private: + int i; + size_t n; + public: + bounded_int(int i, size_t n) : i(i), n(n) {} + + bounded_int& operator--() { + --i; + if (i == -1) { + i = (int)n - 1; + } + return *this; + } + + bounded_int& operator++() { + ++i; + if (i == (int)n) { + i = 0; + } + return *this; + } + + operator int() { return i; } + }; +} + +bool IfcGeom::util::wire_intersections(const TopoDS_Wire& wire, TopTools_ListOfShape& wires, double eps, double eps_real) { + if (!wire.Closed()) { + wires.Append(wire); + return false; + } + + int n = IfcGeom::Kernel::count(wire, TopAbs_EDGE); + if (n < 3) { + wires.Append(wire); + return false; + } + + // Note: initialize empty + Handle(ShapeExtend_WireData) wd = new ShapeExtend_WireData(); + + // ... to be sure to get consecutive edges + BRepTools_WireExplorer exp(wire); + IfcGeom::impl::tree tree; + + int edge_idx = 0; + for (; exp.More(); exp.Next()) { + wd->Add(exp.Current()); + if (n > 64) { + // tfk: indices in tree are 0-based vd 1-based in wiredata + tree.add(edge_idx++, exp.Current()); + } + } + + if (wd->NbEdges() != n) { + // If the number of edges differs, BRepTools_WireExplorer did not + // reach every edge, probably due to loops exactly at vertex locations. + // This is not supported by this algorithm which only elimates loops + // due to edge crossings. + + throw geometry_exception("Invalid loop"); + } + + bool intersected = false; + + // tfk: Extrema on infinite curves proved to be more robust. + // TopoDS_Face face = BRepBuilderAPI_MakeFace(wire, true).Face(); + // ShapeAnalysis_Wire saw(wd, face, getValue(GV_PRECISION)); + + // @todo: should this start from 0 in case of n > 64? + for (int i = 2; i < n; ++i) { + + std::vector js; + if (n > 64) { + Bnd_Box b; + BRepBndLib::Add(wd->Edge(i + 1), b); + b.Enlarge(eps); + js = tree.select_box(b, false); + } else { + boost::push_back(js, boost::irange(0, i - 1)); + } + + for (std::vector::const_iterator it = js.begin(); it != js.end(); ++it) { + int j = *it; + + if (n > 64) { + if (j > i) { + continue; + } + if ((std::max)(i, j) - (std::min)(i, j) <= 1) { + continue; + } + } + + // Only check non-consecutive edges + if (i == n - 1 && j == 0) continue; + + double u11, u12, u21, u22, U1, U2; + GeomAPI_ExtremaCurveCurve ecc( + BRep_Tool::Curve(wd->Edge(i + 1), u11, u12), + BRep_Tool::Curve(wd->Edge(j + 1), u21, u22) + ); + + // @todo: extend this to work in case of multiple extrema and curved segments. + const bool unbounded_intersects = (!ecc.Extrema().IsParallel() && ecc.NbExtrema() == 1 && ecc.Distance(1) < eps); + if (unbounded_intersects) { + ecc.Parameters(1, U1, U2); + + if (u11 > u12) { + std::swap(u11, u12); + } + if (u21 > u22) { + std::swap(u21, u22); + } + + /// @todo: tfk: probably need different thresholds on non-linear curves + u11 -= eps; + u12 += eps; + u21 -= eps; + u22 += eps; + + // tfk: code below is for ShapeAnalysis_Wire::CheckIntersectingEdges() + // IntRes2d_SequenceOfIntersectionPoint points2d; + // TColgp_SequenceOfPnt points3d; + // TColStd_SequenceOfReal errors; + // if (saw.CheckIntersectingEdges(i + 1, j + 1, points2d, points3d, errors)) { + + if (u11 < U1 && U1 < u12 && u21 < U2 && U2 < u22) { + + intersected = true; + + // Explore a forward and backward cycle from the intersection point + for (int fb = 0; fb <= 1; ++fb) { + const bool forward = fb == 0; + + BRepBuilderAPI_MakeWire mw; + bool first = true; + + for (bounded_int k(j, n);;) { + bool intersecting = k == j || k == i; + if (intersecting) { + TopoDS_Edge e = wd->Edge(k + 1); + + TopoDS_Vertex v1, v2; + TopExp::Vertices(e, v1, v2, true); + const TopoDS_Vertex* v = first == forward ? &v2 : &v1; + + // gp_Pnt p2 = points3d.Value(1); + + gp_Pnt p1 = BRep_Tool::Pnt(*v); + gp_Pnt pp1, pp2; + ecc.Points(1, pp1, pp2); + const gp_Pnt& p2 = k == i ? pp1 : pp2; + + // Substitute with a new edge from/to the intersection point + if (p1.Distance(p2) > eps_real * 2) { + double _, __; + Handle_Geom_Curve crv = BRep_Tool::Curve(e, _, __); + BRepBuilderAPI_MakeEdge me(crv, p1, p2); + TopoDS_Edge ed = me.Edge(); + mw.Add(ed); + } + + first = false; + } else { + // Re-use original edge + mw.Add(wd->Edge(k + 1)); + } + + if (k == i) { + break; + } + + if (forward) { + ++k; + } else { + --k; + } + } + + // Recursively process both cuts + + // @todo this is a change in behaviour with eps precomputed from the kernel + // instead of adaptively calculated for the successive iterations. + wire_intersections(mw.Wire(), wires, eps, eps_real); + } + + return true; + } + + } + } + } + + // No intersections found, append original wire + if (!intersected) { + wires.Append(wire); + } + + return intersected; +} + +void IfcGeom::util::select_largest(const TopTools_ListOfShape& shapes, TopoDS_Shape& largest) { + double mass = 0.; + TopTools_ListIteratorOfListOfShape it(shapes); + for (; it.More(); it.Next()) { + /* + // tfk: bounding box is more efficient probably + const TopoDS_Wire& w = TopoDS::Wire(it.Value()); + TopoDS_Face face = BRepBuilderAPI_MakeFace(w).Face(); + const double m = face_area(face); + */ + + Bnd_Box bb; + BRepBndLib::AddClose(it.Value(), bb); + double xyz_min[3], xyz_max[3]; + bb.Get(xyz_min[0], xyz_min[1], xyz_min[2], xyz_max[0], xyz_max[1], xyz_max[2]); + + // @todo hard coded precision. + // @todo this is a really strange measure for wire size. Why not use newell's + // method to project to plane and then calculate size of the 2d bbox? + const double eps = 1.e-5; + + double m = 1.; + for (int i = 0; i < 3; ++i) { + if (Precision::IsNegativeInfinite(xyz_min[i])) { + xyz_min[i] = 0.; + } + if (Precision::IsInfinite(xyz_max[i])) { + xyz_max[i] = 0.; + } + m *= (xyz_max[i] + eps) - (xyz_min[i] - eps); + } + + if (m > mass) { + mass = m; + largest = it.Value(); + } + } +} \ No newline at end of file diff --git a/src/ifcgeom_schema_agnostic/wire_utils.h b/src/ifcgeom_schema_agnostic/wire_utils.h new file mode 100644 index 0000000000..3d905a7338 --- /dev/null +++ b/src/ifcgeom_schema_agnostic/wire_utils.h @@ -0,0 +1,28 @@ +#include +#include +#include + +#include + +namespace IfcGeom { + namespace util { + bool approximate_plane_through_wire(const TopoDS_Wire& wire, gp_Pln& plane, double eps); + + bool flatten_wire(TopoDS_Wire& wire, double eps); + + enum triangulate_wire_result { + TRIANGULATE_WIRE_FAIL, + TRIANGULATE_WIRE_OK, + TRIANGULATE_WIRE_NON_MANIFOLD, + }; + + /// Triangulate the set of wires. The firstmost wire is assumed to be the outer wire. + triangulate_wire_result triangulate_wire(const std::vector& wires, TopTools_ListOfShape& faces); + + // eps: tolerance added to wire intersection checks, can be zero + // eps_real: tolerance used to construct new edge geometry around intersection points, cannot be zero + bool wire_intersections(const TopoDS_Wire& wire, TopTools_ListOfShape& wires, double eps, double eps_real); + + void select_largest(const TopTools_ListOfShape& shapes, TopoDS_Shape& largest); + } +}