From 6f5f18944329f0edf43739d2f5f27c1195bfe44f Mon Sep 17 00:00:00 2001 From: Thomas Krijnen Date: Tue, 14 Jun 2022 15:24:38 +0200 Subject: [PATCH] Missing files --- src/ifcgeom_schema_agnostic/IteratorCache.h | 31 + src/ifcgeom_schema_agnostic/boolean_utils.cpp | 777 ++++++++++++++++++ src/ifcgeom_schema_agnostic/boolean_utils.h | 90 ++ .../face_definition.cpp | 22 + src/ifcgeom_schema_agnostic/face_definition.h | 78 ++ src/ifcgeom_schema_agnostic/sweep_utils.cpp | 356 ++++++++ src/ifcgeom_schema_agnostic/sweep_utils.h | 59 ++ src/ifcgeom_schema_agnostic/wire_builder.cpp | 226 +++++ src/ifcgeom_schema_agnostic/wire_builder.h | 105 +++ 9 files changed, 1744 insertions(+) create mode 100644 src/ifcgeom_schema_agnostic/IteratorCache.h create mode 100644 src/ifcgeom_schema_agnostic/boolean_utils.cpp create mode 100644 src/ifcgeom_schema_agnostic/boolean_utils.h create mode 100644 src/ifcgeom_schema_agnostic/face_definition.cpp create mode 100644 src/ifcgeom_schema_agnostic/face_definition.h create mode 100644 src/ifcgeom_schema_agnostic/sweep_utils.cpp create mode 100644 src/ifcgeom_schema_agnostic/sweep_utils.h create mode 100644 src/ifcgeom_schema_agnostic/wire_builder.cpp create mode 100644 src/ifcgeom_schema_agnostic/wire_builder.h diff --git a/src/ifcgeom_schema_agnostic/IteratorCache.h b/src/ifcgeom_schema_agnostic/IteratorCache.h new file mode 100644 index 0000000000..2f92028359 --- /dev/null +++ b/src/ifcgeom_schema_agnostic/IteratorCache.h @@ -0,0 +1,31 @@ +/******************************************************************************** + * * + * This file is part of IfcOpenShell. * + * * + * IfcOpenShell is free software: you can redistribute it and/or modify * + * it under the terms of the Lesser GNU General Public License as published by * + * the Free Software Foundation, either version 3.0 of the License, or * + * (at your option) any later version. * + * * + * IfcOpenShell is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * Lesser GNU General Public License for more details. * + * * + * You should have received a copy of the Lesser GNU General Public License * + * along with this program. If not, see . * + * * + ********************************************************************************/ + +#ifndef ITERATOR_CACHE_H +#define ITERATOR_CACHE_H + +namespace IfcGeom { + + class IteratorCache { + + }; + +} + +#endif \ No newline at end of file diff --git a/src/ifcgeom_schema_agnostic/boolean_utils.cpp b/src/ifcgeom_schema_agnostic/boolean_utils.cpp new file mode 100644 index 0000000000..3a22488943 --- /dev/null +++ b/src/ifcgeom_schema_agnostic/boolean_utils.cpp @@ -0,0 +1,777 @@ +#include "boolean_utils.h" + +#include "../ifcgeom/IfcGeomTree.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +void IfcGeom::util::copy_operand(const TopTools_ListOfShape & l, TopTools_ListOfShape & r) { +#if OCC_VERSION_HEX < 0x70000 + TopTools_ListIteratorOfListOfShape it(l); + for (; it.More(); it.Next()) { + r.Append(BRepBuilderAPI_Copy(it.Value())); + } +#else + // On OCCT 7.0 and higher BRepAlgoAPI_BuilderAlgo::SetNonDestructive(true) is + // called. Not entirely sure on the behaviour before 7.0, so overcautiously + // create copies. + r.Assign(l); +#endif +} + +TopoDS_Shape IfcGeom::util::copy_operand(const TopoDS_Shape & s) { +#if OCC_VERSION_HEX < 0x70000 + return BRepBuilderAPI_Copy(s); +#else + return s; +#endif +} + +double IfcGeom::util::min_edge_length(const TopoDS_Shape & a) { + double min_edge_len = std::numeric_limits::infinity(); + TopExp_Explorer exp(a, TopAbs_EDGE); + for (; exp.More(); exp.Next()) { + GProp_GProps prop; + BRepGProp::LinearProperties(exp.Current(), prop); + double l = prop.Mass(); + if (l < min_edge_len) { + min_edge_len = l; + } + } + return min_edge_len; +} + +double IfcGeom::util::min_vertex_edge_distance(const TopoDS_Shape & a, double min_search, double max_search) { + double M = std::numeric_limits::infinity(); + + TopTools_IndexedMapOfShape vertices, edges; + + TopExp::MapShapes(a, TopAbs_VERTEX, vertices); + TopExp::MapShapes(a, TopAbs_EDGE, edges); + + IfcGeom::impl::tree tree; + + // Add edges to tree + for (int i = 1; i <= edges.Extent(); ++i) { + tree.add(i, edges(i)); + } + + for (int j = 1; j <= vertices.Extent(); ++j) { + const TopoDS_Vertex& v = TopoDS::Vertex(vertices(j)); + gp_Pnt p = BRep_Tool::Pnt(v); + + Bnd_Box b; + b.Add(p); + b.Enlarge(max_search); + + std::vector edge_idxs = tree.select_box(b, false); + std::vector::const_iterator it = edge_idxs.begin(); + for (; it != edge_idxs.end(); ++it) { + const TopoDS_Edge& e = TopoDS::Edge(edges(*it)); + TopoDS_Vertex v1, v2; + TopExp::Vertices(e, v1, v2); + + if (v.IsSame(v1) || v.IsSame(v2)) { + continue; + } + + BRepAdaptor_Curve crv(e); + Extrema_ExtPC ext(p, crv); + if (!ext.IsDone()) { + continue; + } + + for (int i = 1; i <= ext.NbExt(); ++i) { + const double m = sqrt(ext.SquareDistance(i)); + if (m < M && m > min_search) { + M = m; + } + } + } + } + + return M; +} + +bool IfcGeom::util::faces_overlap(const TopoDS_Face & f, const TopoDS_Face & g) { + points_on_planar_face_generator pgen(f); + + BRep_Builder B; + gp_Pnt test; + double eps = BRep_Tool::Tolerance(f) + BRep_Tool::Tolerance(g); + + BRepExtrema_DistShapeShape x; + x.LoadS1(g); + + while (pgen(test)) { + TopoDS_Vertex V; + B.MakeVertex(V, test, Precision::Confusion()); + x.LoadS2(V); + x.Perform(); + if (x.IsDone() && x.NbSolution() == 1) { + if (x.Value() > eps) { + return false; + } + } + } + + return true; +} + +double IfcGeom::util::min_face_face_distance(const TopoDS_Shape & a, double max_search) { + /* + NB: This is currently only implemented for planar surfaces. + */ + double M = std::numeric_limits::infinity(); + + TopTools_IndexedMapOfShape faces; + + TopExp::MapShapes(a, TopAbs_FACE, faces); + + IfcGeom::impl::tree tree; + + // Add faces to tree + for (int i = 1; i <= faces.Extent(); ++i) { + if (BRep_Tool::Surface(TopoDS::Face(faces(i)))->DynamicType() == STANDARD_TYPE(Geom_Plane)) { + tree.add(i, faces(i)); + } + } + + for (int j = 1; j <= faces.Extent(); ++j) { + const TopoDS_Face& f = TopoDS::Face(faces(j)); + const Handle(Geom_Surface)& fs = BRep_Tool::Surface(f); + + if (fs->DynamicType() != STANDARD_TYPE(Geom_Plane)) { + continue; + } + + points_on_planar_face_generator pgen(f); + + Bnd_Box b; + BRepBndLib::AddClose(f, b); + b.Enlarge(max_search); + + std::vector face_idxs = tree.select_box(b, false); + std::vector::const_iterator it = face_idxs.begin(); + for (; it != face_idxs.end(); ++it) { + if (*it == j) { + continue; + } + + const TopoDS_Face& g = TopoDS::Face(faces(*it)); + const Handle(Geom_Surface)& gs = BRep_Tool::Surface(g); + + auto p0 = Handle(Geom_Plane)::DownCast(fs); + auto p1 = Handle(Geom_Plane)::DownCast(gs); + + if (p0->Position().IsCoplanar(p1->Position(), max_search, asin(max_search))) { + pgen.reset(); + + BRepTopAdaptor_FClass2d cls(g, BRep_Tool::Tolerance(g)); + + gp_Pnt test; + while (pgen(test)) { + gp_Vec d = test.XYZ() - p1->Position().Location().XYZ(); + double u = d.Dot(p1->Position().XDirection()); + double v = d.Dot(p1->Position().YDirection()); + + // nb: TopAbs_ON is explicitly not considered to prevent matching adjacent faces + // with similar orientations. + if (cls.Perform(gp_Pnt2d(u, v)) == TopAbs_IN) { + gp_Pnt test2; + p1->D0(u, v, test2); + double w = std::abs(gp_Vec(p1->Position().Direction().XYZ()).Dot(test2.XYZ() - test.XYZ())); + if (w < M) { + M = w; + } + } + } + } + } + } + + return M; +} + +int IfcGeom::util::bounding_box_overlap(double p, const TopoDS_Shape & a, const TopTools_ListOfShape & b, TopTools_ListOfShape & c) { + int N = 0; + + Bnd_Box A; + BRepBndLib::Add(a, A); + + if (A.IsVoid()) { + return 0; + } + + TopTools_ListIteratorOfListOfShape it(b); + for (; it.More(); it.Next()) { + Bnd_Box B; + BRepBndLib::Add(it.Value(), B); + + if (B.IsVoid()) { + continue; + } + + if (A.Distance(B) < p) { + c.Append(it.Value()); + } else { + ++N; + } + } + + return N; +} + +bool IfcGeom::util::get_edge_axis(const TopoDS_Edge & e, gp_Ax1 & ax) { + double _, __; + + auto crv = BRep_Tool::Curve(e, _, __); + auto line = Handle_Geom_Line::DownCast(crv); + auto bsple = Handle_Geom_BSplineCurve::DownCast(crv); + + if (line) { + ax = line->Position(); + return true; + } else if (bsple) { + if (bsple->NbPoles() == 2 && bsple->Degree() == 1) { + gp_Dir V(bsple->Poles().Last().XYZ() - bsple->Poles().First().XYZ()); + ax = gp_Ax1(bsple->Poles().First(), V); + return true; + } + } + + return false; +} + +bool IfcGeom::util::is_subset(const TopTools_IndexedMapOfShape & lhs, const TopTools_IndexedMapOfShape & rhs) { + if (rhs.Extent() < lhs.Extent()) { + return false; + } + for (int i = 1; i < lhs.Extent(); ++i) { + auto& s = lhs.FindKey(i); + if (!rhs.Contains(s)) { + return false; + } + } + return true; +} + +bool IfcGeom::util::is_extrusion(const gp_Vec & v, const TopoDS_Shape & s, TopoDS_Face & base, std::pair& interval) { + // This assumes UnifySameDomain has been processed on s, so that + // the extrusion top and bottom are a single face. + + TopTools_IndexedDataMapOfShapeListOfShape mapping; + TopExp::MapShapesAndAncestors(s, TopAbs_EDGE, TopAbs_FACE, mapping); + TopExp::MapShapesAndAncestors(s, TopAbs_VERTEX, TopAbs_FACE, mapping); + + TopTools_ListOfShape parallel; + TopTools_IndexedMapOfShape curved_orthogonal; + gp_Ax1 ax; + gp_Ax1 V(gp::Origin(), v); + + // Segment edges in parallel to extrusion direction, and orthogonal or curved, + // where the latter two categories have to make the edges part of the base or + // top face. When neither of these categories the shape is not a extrusion + // or the extrusion direction is not orthogonal to its basis. + for (int i = 1; i < mapping.Extent(); ++i) { + auto& s = mapping.FindKey(i); + if (s.ShapeType() != TopAbs_EDGE) { + continue; + } + + // @todo use a linear tolernace and the face extrimities, see #2218 + const TopoDS_Edge& e = TopoDS::Edge(s); + if (!get_edge_axis(e, ax)) { + // curved + curved_orthogonal.Add(e); + } else if (ax.IsParallel(V, 1.e-7)) { + parallel.Append(e); + } else if (ax.IsNormal(V, 1.e-7)) { + // ortho + curved_orthogonal.Add(e); + } else { + return false; + } + } + + // Select the two faces for which their edges are subsets + // of the ortho/curved edges + TopTools_IndexedMapOfShape ortho_faces; + for (TopExp_Explorer exp(s, TopAbs_FACE); exp.More(); exp.Next()) { + TopTools_IndexedMapOfShape face_edges; + TopExp::MapShapes(exp.Current(), TopAbs_EDGE, face_edges); + if (is_subset(face_edges, curved_orthogonal)) { + ortho_faces.Add(exp.Current()); + } + } + + // There should be a basis and top face + if (ortho_faces.Extent() != 2) { + return false; + } + + // For the parallel edges assert that its two vertices are part + // of both the basis and the top face. + for (TopTools_ListIteratorOfListOfShape it(parallel); + it.More(); it.Next()) { + TopoDS_Vertex v01[2]; + TopExp::Vertices(TopoDS::Edge(it.Value()), v01[0], v01[1]); + + TopTools_IndexedMapOfShape v_ortho_faces; + int nb_ortho_faces[2] = { 0,0 }; + + for (int i = 0; i < 2; ++i) { + auto& faces = mapping.FindFromKey(v01[i]); + + for (TopTools_ListIteratorOfListOfShape jt(faces); + jt.More(); jt.Next()) { + if (ortho_faces.Contains(jt.Value())) { + nb_ortho_faces[i] ++; + v_ortho_faces.Add(jt.Value()); + } + } + } + + bool sets_equal = v_ortho_faces.Size() == ortho_faces.Size() && is_subset(v_ortho_faces, ortho_faces); + if (!sets_equal) { + return false; + } + } + + // Assert the base/top faces are planar and get the interval + // (dot products along axis) for which the extrusion is defined + // If necessary swap the two faces so that the basis face has + // the smallest dot product along the axis. + auto f0 = TopoDS::Face(ortho_faces.FindKey(1)); + auto f1 = TopoDS::Face(ortho_faces.FindKey(2)); + + const Handle(Geom_Surface)& f0_s = BRep_Tool::Surface(f0); + const Handle(Geom_Surface)& f1_s = BRep_Tool::Surface(f1); + + auto p0 = Handle(Geom_Plane)::DownCast(f0_s); + auto p1 = Handle(Geom_Plane)::DownCast(f1_s); + + if (p0.IsNull() || p1.IsNull()) { + return false; + } + + auto dot0 = p0->Location().XYZ().Dot(v.XYZ()); + auto dot1 = p1->Location().XYZ().Dot(v.XYZ()); + + if (dot0 > dot1) { + std::swap(dot0, dot1); + std::swap(f0, f1); + } + + base = f0; + interval = { dot0, dot1 }; + + return true; +} + +int IfcGeom::util::eliminate_touching_operands(double prec, const TopoDS_Shape & a, const TopTools_ListOfShape & bs, TopTools_ListOfShape & c) { + TopTools_IndexedMapOfShape a_faces; + TopExp::MapShapes(a, TopAbs_FACE, a_faces); + + // Check if any of the faces in a are non-planar, which is + // not supported by this quick check. + for (int i = 1; i <= a_faces.Extent(); ++i) { + auto surf = BRep_Tool::Surface(TopoDS::Face(a_faces(i))); + if (surf->DynamicType() != STANDARD_TYPE(Geom_Plane)) { + return 0; + } + } + + TopTools_IndexedMapOfShape a_vertices; + TopExp::MapShapes(a, TopAbs_VERTEX, a_vertices); + + IfcGeom::impl::tree tree; + + // Add faces to tree + for (int i = 1; i <= a_faces.Extent(); ++i) { + tree.add(i, a_faces(i)); + } + + int N = 0; + + TopTools_ListIteratorOfListOfShape it(bs); + for (; it.More(); it.Next()) { + bool is_touching = false; + + auto& b = it.Value(); + + TopTools_IndexedMapOfShape b_faces; + TopExp::MapShapes(b, TopAbs_FACE, b_faces); + + // Check if any of the faces in b are non-planar, which is + // not supported by this quick check. + for (int i = 1; i <= b_faces.Extent(); ++i) { + auto surf = BRep_Tool::Surface(TopoDS::Face(b_faces(i))); + if (surf->DynamicType() != STANDARD_TYPE(Geom_Plane)) { + continue; + } + } + + TopTools_IndexedMapOfShape b_vertices; + TopExp::MapShapes(b, TopAbs_VERTEX, b_vertices); + + for (int k = 1; k <= b_faces.Extent(); ++k) { + const TopoDS_Face& f_b = TopoDS::Face(b_faces(k)); + Bnd_Box B; + BRepBndLib::Add(f_b, B); + + // Query tree using b_face bounding box + for (auto& i : tree.select_box(B, false)) { + const TopoDS_Face& f_a = TopoDS::Face(a_faces(i)); + + TopTools_IndexedMapOfShape f_a_vertices; + TopExp::MapShapes(f_a, TopAbs_VERTEX, f_a_vertices); + + BRepGProp_Face prop_a(f_a); + BRepGProp_Face prop_b(f_b); + + gp_Pnt p_a, p_b; + gp_Vec v_a, v_b; + + double u0, u1, v0, v1; + prop_a.Bounds(u0, u1, v0, v1); + prop_a.Normal((u0 + u1) / 2., (u0 + u1) / 2., p_a, v_a); + + prop_b.Bounds(u0, u1, v0, v1); + prop_b.Normal((u0 + u1) / 2., (u0 + u1) / 2., p_b, v_b); + + bool all_vertices_behind_f_a = true; + + // Check if all 'other' vertices in a are pointing + // away from the face in a, so that there is no geometry + // from a in front of the face that could participate + // in the boolean subtraction. + for (int j = 1; j <= a_vertices.Extent(); ++j) { + if (!f_a_vertices.Contains(a_vertices(j))) { + auto p = BRep_Tool::Pnt(TopoDS::Vertex(a_vertices(j))); + if ((p.XYZ() - p_a.XYZ()).Dot(v_a.XYZ()) > prec) { + all_vertices_behind_f_a = false; + break; + } + } + } + + if (!all_vertices_behind_f_a) { + continue; + } + + // Check if surface normals are opposite + if (v_a.IsOpposite(v_b, 1.e-5)) { + // Check if faces are co-planar + if ((p_b.XYZ() - p_a.XYZ()).Dot(v_a.XYZ()) <= prec) { + + TopTools_IndexedMapOfShape f_b_vertices; + TopExp::MapShapes(f_b, TopAbs_VERTEX, f_b_vertices); + + bool all_vertices_behind_f_b = true; + + // Check if all 'other' vertices in b are pointing + // away from the face in a. So that a boolean subtraction + // would not alter a. + for (int j = 1; j <= b_vertices.Extent(); ++j) { + if (!f_b_vertices.Contains(b_vertices(j))) { + auto p = BRep_Tool::Pnt(TopoDS::Vertex(b_vertices(j))); + if ((p.XYZ() - p_a.XYZ()).Dot(v_a.XYZ()) < prec * 10.) { + all_vertices_behind_f_b = false; + break; + } + } + } + + if (all_vertices_behind_f_b) { + is_touching = true; + break; + } + + } + } + } + + if (is_touching) { + break; + } + } + + if (!is_touching) { + c.Append(it.Value()); + } else { + ++N; + } + } + + return N; +} + +TopoDS_Shape IfcGeom::util::unify(const TopoDS_Shape & s, double tolerance) { + tolerance = (std::min)(min_edge_length(s) / 2., tolerance); + ShapeUpgrade_UnifySameDomain usd(s); +#if OCC_VERSION_HEX >= 0x70200 + usd.SetSafeInputMode(true); +#endif +#if OCC_VERSION_HEX >= 0x70100 + usd.SetLinearTolerance(tolerance); + usd.SetAngularTolerance(1.e-3); +#endif + usd.Build(); + return usd.Shape(); +} + +bool IfcGeom::util::boolean_subtraction_2d_using_builder(const TopoDS_Shape & a_input, const TopTools_ListOfShape & b_input, TopoDS_Shape & result, double eps) { + IfcGeom::impl::tree edge_tree; + + TopTools_ListOfShape ab_input = b_input; + ab_input.Prepend(a_input); + + TopTools_ListIteratorOfListOfShape it(ab_input); + int shape_index = 0; + int edge_index = 0; + std::map edge_index_to_shape_index; + + std::vector shapes; + std::vector> edges; + // First is the outer wire + std::vector wires; + + for (; it.More(); it.Next(), ++shape_index) { + if (it.Value().ShapeType() != TopAbs_FACE) { + return false; + } + + const TopoDS_Face& f = TopoDS::Face(it.Value()); + TopoDS_Wire outer_wire; + + if (shape_index == 0) { + outer_wire = BRepTools::OuterWire(f); + wires.push_back(outer_wire); + } + + size_t num_wires = 0; + TopoDS_Iterator it2(it.Value()); + for (; it2.More(); it2.Next()) { + ++num_wires; + + if (outer_wire.IsNull() || !it2.Value().IsSame(outer_wire)) { + wires.push_back(TopoDS::Wire(it2.Value())); + + if (shape_index == 0 && num_wires > 0) { + // An inner wire on the first operand face: reverse, because + // MakeFace expects inner boundaries to be added as bounded + // areas. + wires.back().Reverse(); + } + } + } + + if (num_wires > 1 && shape_index != 0) { + // The first operand can have inner wires, but the others + // can't because a inner wire would result in an additional + // outer wire for the result. + return false; + } + + shapes.push_back(it.Value()); + TopExp_Explorer exp(it.Value(), TopAbs_EDGE); + for (; exp.More(); exp.Next(), ++edge_index) { + edge_tree.add(edge_index, exp.Current()); + edge_index_to_shape_index[edge_index] = shape_index; + edges.push_back({ shape_index, TopoDS::Edge(exp.Current()) }); + } + } + + shape_index = 0; + edge_index = 0; + + it.Initialize(ab_input); + for (; it.More(); it.Next(), ++shape_index) { + TopExp_Explorer exp(it.Value(), TopAbs_EDGE); + for (; exp.More(); exp.Next(), ++edge_index) { + Bnd_Box b; + BRepBndLib::Add(exp.Current(), b); + b.Enlarge(eps); + + for (auto& i : edge_tree.select_box(b)) { + if (i == edge_index) { + // Skip self-selection + continue; + } + + if (edges[i].first == shape_index) { + // Skip edges of the same operand + continue; + } + + const TopoDS_Edge& e0 = TopoDS::Edge(exp.Current()); + const TopoDS_Edge& e1 = edges[i].second; + + double u11, u12, u21, u22, U1, U2; + + GeomAPI_ExtremaCurveCurve ecc( + BRep_Tool::Curve(e0, u11, u12), + BRep_Tool::Curve(e1, 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; + + if (u11 < U1 && U1 < u12 && u21 < U2 && U2 < u22) { + // Edge curves belonging to different operands intersect, don't process + // using builder. + Logger::Notice("Intersecting boundaries"); + return false; + } + } + } + } + } + + // Only inner wires are considered that are directly contained in the outer wire + // Redundant subtractions are eliminated. + + std::vector redundant(wires.size(), false); + + std::vector wire_faces; + wire_faces.reserve(wires.size()); + + std::vector wire_clss; + wire_clss.reserve(wires.size()); + + std::vector> sass; + sass.reserve(wires.size()); + + for (auto& w : wires) { + wire_faces.push_back(BRepBuilderAPI_MakeFace(w).Face()); + wire_clss.emplace_back(wire_faces.back(), eps); + sass.push_back(std::make_unique(BRep_Tool::Surface(wire_faces.back()))); + } + + // First check for containment in outer wire + for (auto it = ++wires.begin(); it != wires.end(); ++it) { + // Considering a single vertex is sufficient because we have already + // guaranteed that the edges of different operands do not cross. + TopoDS_Iterator it_ed(*it); + auto& ed = it_ed.Value(); + + TopoDS_Iterator it_v(ed); + auto& v = TopoDS::Vertex(it_v.Value()); + + auto pnt = BRep_Tool::Pnt(v); + auto p2d = sass[0]->ValueOfUV(pnt, eps); + if (wire_clss[0].Perform(p2d) != TopAbs_IN) { + // A wire is not contained in the outer wire, it's a subtraction without + // any effect and marked as redundant. Feeding it to the builder algo + // will likely cause problems. + redundant[std::distance(wires.begin(), it)] = true; + Logger::Notice("Subtraction operand outside of outer bound"); + } + } + + // Now build a tree to find inner wires contained in other inner wires + // NB first wire is *not* in this tree + IfcGeom::impl::tree wire_tree; + for (size_t wire_index = 1; wire_index < wires.size(); ++wire_index) { + wire_tree.add(wire_index, wires[wire_index]); + } + + for (size_t wire_index = 1; wire_index < wires.size(); ++wire_index) { + Bnd_Box b; + BRepBndLib::Add(wires[wire_index], b); + b.Enlarge(eps); + + // We're only selecting operands completely within b because we + // have already guaranteed they do not intersect. So they are + // either fully in or out. Selecting with complete_within=true + // will filter out some unnecessary cases. It also means we need + // that due this assymetry we need to process all pairs of wire + // indices and not just the pairs where the first element is less + // than the second element. + for (auto& other_index : wire_tree.select_box(b, true)) { + // other_index is fully contained in wire_index + if (wire_index == other_index) { + continue; + } + + TopoDS_Iterator it_ed(wires[other_index]); + auto& ed = it_ed.Value(); + + TopoDS_Iterator it_v(ed); + auto& v = TopoDS::Vertex(it_v.Value()); + + auto pnt = BRep_Tool::Pnt(v); + auto p2d = sass[wire_index]->ValueOfUV(pnt, eps); + if (wire_clss[wire_index].Perform(p2d) == TopAbs_IN) { + // A wire is contained within another operand + redundant[other_index] = true; + Logger::Notice("Subtraction operand contained in other"); + } + } + } + + BRepBuilderAPI_MakeFace mf(wire_faces[0]); + for (size_t wire_index = 1; wire_index < wires.size(); ++wire_index) { + if (!redundant[wire_index]) { + mf.Add(TopoDS::Wire(wires[wire_index].Reversed())); + } + } + result = mf.Face(); + + return true; +} + +void IfcGeom::util::points_on_planar_face_generator::reset() { + i = j = (int)inset_; +} + +bool IfcGeom::util::points_on_planar_face_generator::operator()(gp_Pnt& p) { + while (j < N) { + double u = u0 + (u1 - u0) * i / N; + double v = v0 + (v1 - v0) * j / N; + + i++; + if (i == N) { + i = 0; + j++; + } + + // Specifically does not consider ON + if (cls_.Perform(gp_Pnt2d(u, v)) == TopAbs_IN) { + plane_->D0(u, v, p); + return true; + } + } + + return false; +} \ No newline at end of file diff --git a/src/ifcgeom_schema_agnostic/boolean_utils.h b/src/ifcgeom_schema_agnostic/boolean_utils.h new file mode 100644 index 0000000000..31cc9ed4f2 --- /dev/null +++ b/src/ifcgeom_schema_agnostic/boolean_utils.h @@ -0,0 +1,90 @@ +/******************************************************************************** + * * + * This file is part of IfcOpenShell. * + * * + * IfcOpenShell is free software: you can redistribute it and/or modify * + * it under the terms of the Lesser GNU General Public License as published by * + * the Free Software Foundation, either version 3.0 of the License, or * + * (at your option) any later version. * + * * + * IfcOpenShell is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * Lesser GNU General Public License for more details. * + * * + * You should have received a copy of the Lesser GNU General Public License * + * along with this program. If not, see . * + * * + ********************************************************************************/ + +#ifndef BOOLEAN_UTILS_H +#define BOOLEAN_UTILS_H + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace IfcGeom { + namespace util { + + void copy_operand(const TopTools_ListOfShape& l, TopTools_ListOfShape& r); + + TopoDS_Shape copy_operand(const TopoDS_Shape& s); + + double min_edge_length(const TopoDS_Shape& a); + + double min_vertex_edge_distance(const TopoDS_Shape& a, double min_search, double max_search); + + class points_on_planar_face_generator { + private: + const TopoDS_Face& f_; + Handle(Geom_Surface) plane_; + BRepTopAdaptor_FClass2d cls_; + double u0, u1, v0, v1; + int i, j; + bool inset_; + static const int N = 10; + + public: + points_on_planar_face_generator(const TopoDS_Face& f, bool inset = false) + : f_(f) + , plane_(BRep_Tool::Surface(f_)) + , cls_(f_, BRep_Tool::Tolerance(f_)) + , i((int)inset), j((int)inset) + , inset_(inset) + { + BRepTools::UVBounds(f_, u0, u1, v0, v1); + } + + void reset(); + + bool operator()(gp_Pnt& p); + }; + + bool faces_overlap(const TopoDS_Face& f, const TopoDS_Face& g); + + double min_face_face_distance(const TopoDS_Shape& a, double max_search); + + int bounding_box_overlap(double p, const TopoDS_Shape& a, const TopTools_ListOfShape& b, TopTools_ListOfShape& c); + + bool get_edge_axis(const TopoDS_Edge& e, gp_Ax1& ax); + + bool is_subset(const TopTools_IndexedMapOfShape& lhs, const TopTools_IndexedMapOfShape& rhs); + + bool is_extrusion(const gp_Vec& v, const TopoDS_Shape& s, TopoDS_Face& base, std::pair& interval); + + int eliminate_touching_operands(double prec, const TopoDS_Shape& a, const TopTools_ListOfShape& bs, TopTools_ListOfShape& c); + + TopoDS_Shape unify(const TopoDS_Shape& s, double tolerance); + + bool boolean_subtraction_2d_using_builder(const TopoDS_Shape& a_input, const TopTools_ListOfShape& b_input, TopoDS_Shape& result, double eps); + + } +} + +#endif \ No newline at end of file diff --git a/src/ifcgeom_schema_agnostic/face_definition.cpp b/src/ifcgeom_schema_agnostic/face_definition.cpp new file mode 100644 index 0000000000..4111c6c675 --- /dev/null +++ b/src/ifcgeom_schema_agnostic/face_definition.cpp @@ -0,0 +1,22 @@ +#include "face_definition.h" + +#include +#include +#include +#include + +/* Returns whether wire conforms to a polyhedron, i.e. only edges with linear curves*/ +bool IfcGeom::util::is_polyhedron(const TopoDS_Wire & wire) { + double a, b; + TopLoc_Location l; + + TopoDS_Iterator it(wire, false, false); + for (; it.More(); it.Next()) { + auto crv = BRep_Tool::Curve(TopoDS::Edge(it.Value()), l, a, b); + if (!crv || crv->DynamicType() != STANDARD_TYPE(Geom_Line)) { + return false; + } + } + + return true; +} diff --git a/src/ifcgeom_schema_agnostic/face_definition.h b/src/ifcgeom_schema_agnostic/face_definition.h new file mode 100644 index 0000000000..6fe982055b --- /dev/null +++ b/src/ifcgeom_schema_agnostic/face_definition.h @@ -0,0 +1,78 @@ +/******************************************************************************** + * * + * This file is part of IfcOpenShell. * + * * + * IfcOpenShell is free software: you can redistribute it and/or modify * + * it under the terms of the Lesser GNU General Public License as published by * + * the Free Software Foundation, either version 3.0 of the License, or * + * (at your option) any later version. * + * * + * IfcOpenShell is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * Lesser GNU General Public License for more details. * + * * + * You should have received a copy of the Lesser GNU General Public License * + * along with this program. If not, see . * + * * + ********************************************************************************/ + +#ifndef FACE_DEFINITION_H +#define FACE_DEFINITION_H + +#include +#include + +#include +#include + +namespace IfcGeom { + namespace util { + + /* Returns whether wire conforms to a polyhedron, i.e. only edges with linear curves*/ + bool is_polyhedron(const TopoDS_Wire& wire); + + /* A temporary structure to store the intermediate data for the face conversion */ + class face_definition { + private: + Handle(Geom_Surface) surface_; + std::vector wires_; + bool all_outer_; + public: + face_definition() : surface_(), all_outer_(false) {} + + typedef std::vector::const_iterator wire_it; + + bool& all_outer() { + return all_outer_; + } + + bool all_outer() const { + return all_outer_; + } + + Handle(Geom_Surface)& surface() { + return surface_; + } + + const Handle(Geom_Surface)& surface() const { + return surface_; + } + + std::vector& wires() { + return wires_; + } + + const TopoDS_Wire& outer_wire() const { + return wires_.front(); + } + + std::pair inner_wires() const { + return { wires_.begin() + 1, wires_.end() }; + } + }; + + } +} + +#endif diff --git a/src/ifcgeom_schema_agnostic/sweep_utils.cpp b/src/ifcgeom_schema_agnostic/sweep_utils.cpp new file mode 100644 index 0000000000..8c166b7dbb --- /dev/null +++ b/src/ifcgeom_schema_agnostic/sweep_utils.cpp @@ -0,0 +1,356 @@ +#include "sweep_utils.h" + +#include "../ifcparse/IfcLogger.h" +#include "../ifcgeom_schema_agnostic/Kernel.h" + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +bool IfcGeom::util::wire_is_c1_continuous(const TopoDS_Wire & w, double tol) { + // NB Note that c0 continuity is NOT checked! + + TopTools_IndexedDataMapOfShapeListOfShape map; + TopExp::MapShapesAndAncestors(w, TopAbs_VERTEX, TopAbs_EDGE, map); + for (int i = 1; i <= map.Extent(); ++i) { + const auto& li = map.FindFromIndex(i); + if (li.Extent() == 2) { + const TopoDS_Vertex& v = TopoDS::Vertex(map.FindKey(i)); + + const TopoDS_Edge& e0 = TopoDS::Edge(li.First()); + const TopoDS_Edge& e1 = TopoDS::Edge(li.Last()); + + double u0 = BRep_Tool::Parameter(v, e0); + double u1 = BRep_Tool::Parameter(v, e1); + + double _, __; + Handle(Geom_Curve) c0 = BRep_Tool::Curve(e0, _, __); + Handle(Geom_Curve) c1 = BRep_Tool::Curve(e1, _, __); + + gp_Pnt p; + gp_Vec v0, v1; + c0->D1(u0, p, v0); + c1->D1(u1, p, v1); + + if (1. - std::abs(v0.Normalized().Dot(v1.Normalized())) > tol) { + return false; + } + } + } + return true; +} + +bool IfcGeom::util::wire_to_ax(const TopoDS_Wire & wire, gp_Ax2 & directrix) { + gp_Pnt directrix_origin; + gp_Vec directrix_tangent; + + TopoDS_Edge edge; + + // Find first edge + TopoDS_Vertex v0, v1; + TopExp::Vertices(wire, v0, v1); + TopTools_IndexedDataMapOfShapeListOfShape map; + TopExp::MapShapesAndAncestors(wire, TopAbs_VERTEX, TopAbs_EDGE, map); + if (v0.IsSame(v1) && map.Contains(v0) && map.FindFromKey(v0).Extent() == 2) { + // Closed wire, with more than 1 edges + auto es = map.FindFromKey(v0); + auto e1 = TopoDS::Edge(es.First()); + auto e2 = TopoDS::Edge(es.Last()); + + double u0, u1; + + gp_Vec accum; + + Handle(Geom_Curve) crv = BRep_Tool::Curve(e1, u0, u1); + crv->D1(TopExp::FirstVertex(e1).IsSame(v0) ? u0 : u1, directrix_origin, directrix_tangent); + + accum += directrix_tangent; + + crv = BRep_Tool::Curve(e2, u0, u1); + crv->D1(TopExp::FirstVertex(e2).IsSame(v0) ? u0 : u1, directrix_origin, directrix_tangent); + + accum += directrix_tangent; + + directrix_tangent = accum; + + } else if (map.Contains(v0) && map.FindFromKey(v0).Extent() == 1) { + edge = TopoDS::Edge(map.FindFromKey(v0).First()); + + double u0, u1; + Handle(Geom_Curve) crv = BRep_Tool::Curve(edge, u0, u1); + crv->D1(u0, directrix_origin, directrix_tangent); + } else { + Logger::Error("Unable to locate first edge"); + return false; + } + + directrix = gp_Ax2(directrix_origin, directrix_tangent); + + return true; +} + +bool IfcGeom::util::is_single_linear_edge(const TopoDS_Wire & wire) { + TopExp_Explorer exp(wire, TopAbs_EDGE); + if (!exp.More()) { + return false; + } + TopoDS_Edge e = TopoDS::Edge(exp.Current()); + exp.Next(); + if (exp.More()) { + return false; + } + double u, v; + Handle_Geom_Curve crv = BRep_Tool::Curve(e, u, v); + return crv->DynamicType() == STANDARD_TYPE(Geom_Line); +} + +bool IfcGeom::util::is_single_circular_edge(const TopoDS_Wire & wire) { + TopExp_Explorer exp(wire, TopAbs_EDGE); + if (!exp.More()) { + return false; + } + TopoDS_Edge e = TopoDS::Edge(exp.Current()); + exp.Next(); + if (exp.More()) { + return false; + } + double u, v; + Handle_Geom_Curve crv = BRep_Tool::Curve(e, u, v); + return crv->DynamicType() == STANDARD_TYPE(Geom_Circle); +} + +void IfcGeom::util::process_sweep_as_extrusion(const TopoDS_Wire & wire, const TopoDS_Wire & section, TopoDS_Shape & result) { + TopExp_Explorer exp(wire, TopAbs_EDGE); + TopoDS_Edge e = TopoDS::Edge(exp.Current()); + double u, v; + Handle_Geom_Curve crv = BRep_Tool::Curve(e, u, v); + const auto& dir = Handle(Geom_Line)::DownCast(crv)->Position().Direction(); + // OCCT line is normalized so diff in parametric coords equals length + const double depth = std::abs(u - v); + // @todo we could be extruding the wire only when we know this is an intermediate edge. + TopoDS_Face face = BRepBuilderAPI_MakeFace(section).Face(); + result = BRepPrimAPI_MakePrism(face, depth*dir).Shape(); +} + +void IfcGeom::util::process_sweep_as_revolution(const TopoDS_Wire & wire, const TopoDS_Wire & section, TopoDS_Shape & result) { + TopExp_Explorer exp(wire, TopAbs_EDGE); + TopoDS_Edge e = TopoDS::Edge(exp.Current()); + double u, v; + Handle_Geom_Curve crv = BRep_Tool::Curve(e, u, v); + auto circ = Handle(Geom_Circle)::DownCast(crv); + // @todo we could be extruding the wire only when we know this is an intermediate edge. + const double depth = std::abs(u - v); + TopoDS_Face face = BRepBuilderAPI_MakeFace(section).Face(); + result = BRepPrimAPI_MakeRevol(face, circ->Axis(), depth).Shape(); +} + +void IfcGeom::util::process_sweep_as_pipe(const TopoDS_Wire & wire, const TopoDS_Wire & section, TopoDS_Shape & result, bool force_transformed) { + // This tolerance is fairly high due to the linear edge substitution for small (or large radii) conical curves. + const bool is_continuous = wire_is_c1_continuous(wire, 1.e-2); + BRepOffsetAPI_MakePipeShell builder(wire); + builder.Add(section); + builder.SetTransitionMode(is_continuous || force_transformed ? BRepBuilderAPI_Transformed : BRepBuilderAPI_RightCorner); + try { + builder.Build(); + } catch (Standard_Failure& e) { + // We fallback to BRepBuilderAPI_Transformed, but likely with visual artefacts. + if (!(is_continuous || force_transformed)) { + return process_sweep_as_pipe(wire, section, result, true); + } else { + throw e; + } + } + builder.MakeSolid(); + result = builder.Shape(); +} + +void IfcGeom::util::sort_edges(const TopoDS_Wire & wire, std::vector& sorted_edges) { + TopTools_IndexedDataMapOfShapeListOfShape map; + TopExp::MapShapesAndAncestors(wire, TopAbs_VERTEX, TopAbs_EDGE, map); + + for (int i = 1; i <= map.Extent(); ++i) { + if (map.FindFromIndex(i).Extent() > 2) { + Logger::Warning("Self-intersecting Directrix"); + } + } + + std::set seen; + + auto num_edges = IfcGeom::Kernel::count(wire, TopAbs_EDGE); + + TopoDS_Vertex v0, v1; + // @todo this creates the ancestor map twice + TopExp::Vertices(wire, v0, v1); + + bool ignore_first_equality_because_closed = v0.IsSame(v1); + + // @todo this probably still does not work on a closed wire consisting of one (circular) edge. + + while ((int)sorted_edges.size() < num_edges && + (!v0.IsSame(v1) || ignore_first_equality_because_closed)) { + ignore_first_equality_because_closed = false; + if (!map.Contains(v0)) { + throw std::runtime_error("Disconnected vertex"); + } + const TopTools_ListOfShape& es = map.FindFromKey(v0); + TopoDS_Vertex ve0, ve1; + TopTools_ListIteratorOfListOfShape it(es); + bool added = false; + for (; it.More(); it.Next()) { + const TopoDS_Edge& e = TopoDS::Edge(it.Value()); + TopExp::Vertices(e, ve0, ve1, true); + if (ve0.IsSame(v0) && seen.find(&*e.TShape()) == seen.end()) { + sorted_edges.push_back(e); + v0 = ve1; + added = true; + seen.insert(&*e.TShape()); + break; + } + } + if (!added) { + throw std::runtime_error("Disconnected edge"); + } + } +} + +// #939: a closed loop causes failed triangulation in 7.3 and artefacts +// in 7.4 so we break up a closed wire into two equal parts. +void IfcGeom::util::break_closed(const TopoDS_Wire & wire, std::vector& wires) { + std::vector sorted_edges; + sort_edges(wire, sorted_edges); + + if (sorted_edges.size() == 1) { + wires.push_back(wire); + return; + } + + BRep_Builder B; + + wires.emplace_back(); + B.MakeWire(wires.back()); + + for (size_t i = 0; i < sorted_edges.size(); ++i) { + if (i == sorted_edges.size() / 2) { + wires.emplace_back(); + B.MakeWire(wires.back()); + } + + const auto& e = sorted_edges[i]; + B.Add(wires.back(), e); + } +} + +void IfcGeom::util::segment_adjacent_non_linear(const TopoDS_Wire & wire, std::vector& wires) { + std::vector sorted_edges; + sort_edges(wire, sorted_edges); + + BRep_Builder B; + double u, v; + + wires.emplace_back(); + B.MakeWire(wires.back()); + + for (int i = 0; i < (int)sorted_edges.size() - 1; ++i) { + const auto& e = sorted_edges[i]; + Handle_Geom_Curve crv = BRep_Tool::Curve(e, u, v); + const bool is_linear = crv->DynamicType() == STANDARD_TYPE(Geom_Line); + + const auto& f = sorted_edges[i + 1]; + crv = BRep_Tool::Curve(f, u, v); + const bool next_is_linear = crv->DynamicType() == STANDARD_TYPE(Geom_Line); + + B.Add(wires.back(), e); + + if (!is_linear && !next_is_linear) { + wires.emplace_back(); + B.MakeWire(wires.back()); + } + } + + if (!sorted_edges.empty()) { + B.Add(wires.back(), sorted_edges.back()); + } +} + +// @todo make this generic for other sweeps not just swept disk +void IfcGeom::util::process_sweep(const TopoDS_Wire & wire, double radius, TopoDS_Shape & result) { + std::vector wires, wires_tmp; + segment_adjacent_non_linear(wire, wires_tmp); + for (auto& w : wires_tmp) { + break_closed(w, wires); + } + + TopoDS_Compound C; + BRep_Builder B; + if (wires.size() > 1) { + B.MakeCompound(C); + } + + for (auto& w : wires) { + TopoDS_Shape part; + + gp_Ax2 directrix; + if (!wire_to_ax(w, directrix)) { + continue; + } + Handle(Geom_Circle) circle = new Geom_Circle(directrix, radius); + TopoDS_Wire section = BRepBuilderAPI_MakeWire(BRepBuilderAPI_MakeEdge(circle)); + + if (is_single_circular_edge(w)) { + process_sweep_as_revolution(w, section, part); + } else if (is_single_linear_edge(w)) { + process_sweep_as_extrusion(w, section, part); + } else { + process_sweep_as_pipe(w, section, part); + } + if (wires.size() > 1) { + B.Add(C, part); + } else { + result = part; + } + } + + if (wires.size() > 1) { + result = C; + } + + /* + // Eliminate Swept Surfaces? + result = ShapeCustom::SweptToElementary(result); + + // Eliminate Trimmed Surfaces? + ShapeBuild_ReShape sbrs; + BRep_Builder b; + TopExp_Explorer exp(result, TopAbs_FACE); + for (; exp.More(); exp.Next()) { + const TopoDS_Face& f = TopoDS::Face(exp.Current()); + auto S = BRep_Tool::Surface(f); + if (S->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) { + auto RTS = Handle(Geom_RectangularTrimmedSurface)::DownCast(S); + auto B = RTS->BasisSurface(); + TopoDS_Shape newf = f.EmptyCopied(); + // @todo Is it ok to assume no location? + b.MakeFace(TopoDS::Face(newf), B, BRep_Tool::Tolerance(f)); + sbrs.Replace(f, newf); + } + } + result = sbrs.Apply(result); + */ +} diff --git a/src/ifcgeom_schema_agnostic/sweep_utils.h b/src/ifcgeom_schema_agnostic/sweep_utils.h new file mode 100644 index 0000000000..033891d9c9 --- /dev/null +++ b/src/ifcgeom_schema_agnostic/sweep_utils.h @@ -0,0 +1,59 @@ +/******************************************************************************** + * * + * This file is part of IfcOpenShell. * + * * + * IfcOpenShell is free software: you can redistribute it and/or modify * + * it under the terms of the Lesser GNU General Public License as published by * + * the Free Software Foundation, either version 3.0 of the License, or * + * (at your option) any later version. * + * * + * IfcOpenShell is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * Lesser GNU General Public License for more details. * + * * + * You should have received a copy of the Lesser GNU General Public License * + * along with this program. If not, see . * + * * + ********************************************************************************/ + +#ifndef SWEEP_UTILS_H +#define SWEEP_UTILS_H + +#include +#include + +#include + +namespace IfcGeom { + namespace util { + + bool wire_is_c1_continuous(const TopoDS_Wire& w, double tol); + + bool wire_to_ax(const TopoDS_Wire& wire, gp_Ax2& directrix); + + bool is_single_linear_edge(const TopoDS_Wire& wire); + + bool is_single_circular_edge(const TopoDS_Wire& wire); + + void process_sweep_as_extrusion(const TopoDS_Wire& wire, const TopoDS_Wire& section, TopoDS_Shape& result); + + void process_sweep_as_revolution(const TopoDS_Wire& wire, const TopoDS_Wire& section, TopoDS_Shape& result); + + void process_sweep_as_pipe(const TopoDS_Wire& wire, const TopoDS_Wire& section, TopoDS_Shape& result, bool force_transformed = false); + + void sort_edges(const TopoDS_Wire& wire, std::vector& sorted_edges); + + + // #939: a closed loop causes failed triangulation in 7.3 and artefacts + // in 7.4 so we break up a closed wire into two equal parts. + void break_closed(const TopoDS_Wire& wire, std::vector& wires); + + void segment_adjacent_non_linear(const TopoDS_Wire& wire, std::vector& wires); + + // @todo make this generic for other sweeps not just swept disk + void process_sweep(const TopoDS_Wire& wire, double radius, TopoDS_Shape& result); + } +} + +#endif diff --git a/src/ifcgeom_schema_agnostic/wire_builder.cpp b/src/ifcgeom_schema_agnostic/wire_builder.cpp new file mode 100644 index 0000000000..f63f88f498 --- /dev/null +++ b/src/ifcgeom_schema_agnostic/wire_builder.cpp @@ -0,0 +1,226 @@ +#include "wire_builder.h" + +#include "../ifcparse/IfcLogger.h" +#include "../ifcgeom_schema_agnostic/Kernel.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Returns the first edge of a wire +TopoDS_Edge IfcGeom::util::first_edge(const TopoDS_Wire & w) { + TopoDS_Vertex v1, v2; + TopExp::Vertices(w, v1, v2); + TopTools_IndexedDataMapOfShapeListOfShape wm; + TopExp::MapShapesAndAncestors(w, TopAbs_VERTEX, TopAbs_EDGE, wm); + return TopoDS::Edge(wm.FindFromKey(v1).First()); +} + +// Returns new wire with the edge replaced by a linear edge with the vertex v moved to p +TopoDS_Wire IfcGeom::util::adjust(const TopoDS_Wire & w, const TopoDS_Vertex & v, const gp_Pnt & p) { + TopTools_IndexedDataMapOfShapeListOfShape map; + TopExp::MapShapesAndAncestors(w, TopAbs_VERTEX, TopAbs_EDGE, map); + + bool all_linear = true, single_circle = false, first = true; + + const TopTools_ListOfShape& edges = map.FindFromKey(v); + TopTools_ListIteratorOfListOfShape it(edges); + for (; it.More(); it.Next()) { + const TopoDS_Edge& e = TopoDS::Edge(it.Value()); + double _, __; + Handle(Geom_Curve) crv = BRep_Tool::Curve(e, _, __); + const bool is_line = crv->DynamicType() == STANDARD_TYPE(Geom_Line); + const bool is_circle = crv->DynamicType() == STANDARD_TYPE(Geom_Circle); + all_linear = all_linear && is_line; + single_circle = first && is_circle; + } + + if (all_linear) { + BRep_Builder b; + TopoDS_Vertex v2; + b.MakeVertex(v2, p, BRep_Tool::Tolerance(v)); + + ShapeBuild_ReShape reshape; + reshape.Replace(v.Oriented(TopAbs_FORWARD), v2); + + return TopoDS::Wire(reshape.Apply(w)); + } else if (single_circle) { + TopoDS_Vertex v1, v2; + TopExp::Vertices(w, v1, v2); + + gp_Pnt p1, p2, p3; + p1 = v.IsEqual(v1) ? p : BRep_Tool::Pnt(v1); + p3 = v.IsEqual(v2) ? p : BRep_Tool::Pnt(v2); + + double a, b; + Handle(Geom_Curve) crv = BRep_Tool::Curve(TopoDS::Edge(edges.First()), a, b); + crv->D0((a + b) / 2., p2); + + GC_MakeCircle mc(p1, p2, p3); + if (!mc.IsDone()) { + throw IfcGeom::geometry_exception("Failed to adjust circle"); + } + + TopoDS_Edge edge = BRepBuilderAPI_MakeEdge(mc.Value(), p1, p3).Edge(); + BRepBuilderAPI_MakeWire builder; + builder.Add(edge); + return builder.Wire(); + } else { + throw IfcGeom::geometry_exception("Unexpected wire to adjust"); + } +} + +double IfcGeom::util::deflection_for_approximating_circle(double radius, double param) { + return -radius * std::cos(1. / 2. * param) * std::cos(param) - radius * std::sin(1. / 2. * param) * std::sin(param) + radius; +} + +bool IfcGeom::util::create_edge_over_curve_with_log_messages(const Handle_Geom_Curve & crv, const double eps, const gp_Pnt & p1, const gp_Pnt & p2, TopoDS_Edge & result) { + if (crv->IsClosed() && p1.Distance(p2) <= eps) { + BRepBuilderAPI_MakeEdge me(crv); + if (me.IsDone()) { + result = me.Edge(); + return true; + } else { + return false; + } + } + + BRep_Builder builder; + TopoDS_Vertex v1, v2; + /// @todo project first and emit warnings accordingly + builder.MakeVertex(v1, p1, eps); + builder.MakeVertex(v2, p2, eps); + + BRepBuilderAPI_MakeEdge me(crv, v1, v2); + if (!me.IsDone()) { + const double eps2 = eps * eps; + if (me.Error() == BRepBuilderAPI_PointProjectionFailed) { + GeomAdaptor_Curve GAC(crv); + const gp_Pnt* ps[2] = { &p1, &p2 }; + for (int i = 0; i < 2; ++i) { + Extrema_ExtPC extrema(*ps[i], GAC); + if (extrema.IsDone()) { + int n = extrema.NbExt(); + double dmin = std::numeric_limits::infinity(); + for (int j = 1; j <= n; j++) { + const double d = extrema.SquareDistance(j); + if (d < dmin) { + dmin = d; + } + } + if (dmin == std::numeric_limits::infinity()) { + Logger::Error("No extrema for point"); + } else if (dmin > eps2) { + Logger::Error("Distance of " + boost::lexical_cast(std::sqrt(dmin)) + " exceeds tolerance"); + } + } else { + Logger::Error("Failed to calculate extrema for point"); + } + } + } + return false; + } + result = me.Edge(); + return true; +} + +void IfcGeom::util::wire_builder::operator()(const TopoDS_Shape& a) { + const TopoDS_Wire& w = TopoDS::Wire(a); + if (override_next_) { + override_next_ = false; + TopoDS_Edge e = first_edge(w); + mw_.Add(adjust(w, TopExp::FirstVertex(e, true), next_override_)); + } else { + mw_.Add(w); + } +} + +void IfcGeom::util::wire_builder::operator()(const TopoDS_Shape& a, const TopoDS_Shape& b, bool last) { + TopoDS_Wire w1 = TopoDS::Wire(a); + const TopoDS_Wire& w2 = TopoDS::Wire(b); + + if (override_next_) { + override_next_ = false; + TopoDS_Edge e = first_edge(w1); + w1 = adjust(w1, TopExp::FirstVertex(e, true), next_override_); + } + + TopoDS_Vertex w11, w12, w21, w22; + TopExp::Vertices(w1, w11, w12); + TopExp::Vertices(w2, w21, w22); + + gp_Pnt p1 = BRep_Tool::Pnt(w12); + gp_Pnt p2 = BRep_Tool::Pnt(w21); + + double dist = p1.Distance(p2); + + // Distance is within tolerance, this is fine + if (dist < p_) { + mw_.Add(w1); + goto check; + } + + // Distance is too large for attempting to move end points, add intermediate edge + if (dist > 1000. * p_) { + mw_.Add(w1); + mw_.Add(BRepBuilderAPI_MakeEdge(p1, p2)); + Logger::Warning("Added additional segment to close gap with length " + boost::lexical_cast(dist) + " to:", inst_); + goto check; + } + + { + TopTools_IndexedDataMapOfShapeListOfShape wmap1, wmap2; + + // Find edges connected to end- and begin vertex + TopExp::MapShapesAndAncestors(w1, TopAbs_VERTEX, TopAbs_EDGE, wmap1); + TopExp::MapShapesAndAncestors(w2, TopAbs_VERTEX, TopAbs_EDGE, wmap2); + + const TopTools_ListOfShape& last_edges = wmap1.FindFromKey(w12); + const TopTools_ListOfShape& first_edges = wmap2.FindFromKey(w21); + + double _, __; + if (last_edges.Extent() == 1 && first_edges.Extent() == 1) { + Handle(Geom_Curve) c1 = BRep_Tool::Curve(TopoDS::Edge(last_edges.First()), _, __); + Handle(Geom_Curve) c2 = BRep_Tool::Curve(TopoDS::Edge(first_edges.First()), _, __); + + const bool is_line1 = c1->DynamicType() == STANDARD_TYPE(Geom_Line); + const bool is_line2 = c2->DynamicType() == STANDARD_TYPE(Geom_Line); + + const bool is_circle1 = c1->DynamicType() == STANDARD_TYPE(Geom_Circle); + const bool is_circle2 = c2->DynamicType() == STANDARD_TYPE(Geom_Circle); + + // Preferably adjust the segment that is linear + if (is_line1 || (is_circle1 && !is_line2)) { + mw_.Add(adjust(w1, w12, p2)); + Logger::Notice("Adjusted edge end-point with distance " + boost::lexical_cast(dist) + " on:", inst_); + } else if ((is_line2 || is_circle2) && !last) { + mw_.Add(w1); + override_next_ = true; + next_override_ = p1; + Logger::Notice("Adjusted edge end-point with distance " + boost::lexical_cast(dist) + " on:", inst_); + } else { + // In all other cases an edge is added + mw_.Add(w1); + mw_.Add(BRepBuilderAPI_MakeEdge(p1, p2)); + Logger::Warning("Added additional segment to close gap with length " + boost::lexical_cast(dist) + " to:", inst_); + } + } else { + Logger::Error("Internal error, inconsistent wire segments", inst_); + mw_.Add(w1); + } + } + +check: + if (mw_.Error() == BRepBuilderAPI_NonManifoldWire) { + Logger::Error("Non-manifold curve segments:", inst_); + } else if (mw_.Error() == BRepBuilderAPI_DisconnectedWire) { + Logger::Error("Failed to join curve segments:", inst_); + } +} \ No newline at end of file diff --git a/src/ifcgeom_schema_agnostic/wire_builder.h b/src/ifcgeom_schema_agnostic/wire_builder.h new file mode 100644 index 0000000000..7a75c32a7a --- /dev/null +++ b/src/ifcgeom_schema_agnostic/wire_builder.h @@ -0,0 +1,105 @@ +/******************************************************************************** + * * + * This file is part of IfcOpenShell. * + * * + * IfcOpenShell is free software: you can redistribute it and/or modify * + * it under the terms of the Lesser GNU General Public License as published by * + * the Free Software Foundation, either version 3.0 of the License, or * + * (at your option) any later version. * + * * + * IfcOpenShell is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * Lesser GNU General Public License for more details. * + * * + * You should have received a copy of the Lesser GNU General Public License * + * along with this program. If not, see . * + * * + ********************************************************************************/ + +#ifndef WIRE_BUILDER_H +#define WIRE_BUILDER_H + +#include "../ifcparse/IfcBaseClass.h" + +#include + +#include +#include +#include + +#include +#include + +namespace IfcGeom { + namespace util { + // Returns the first edge of a wire + TopoDS_Edge first_edge(const TopoDS_Wire& w); + + // Returns new wire with the edge replaced by a linear edge with the vertex v moved to p + TopoDS_Wire adjust(const TopoDS_Wire& w, const TopoDS_Vertex& v, const gp_Pnt& p); + + // A wrapper around BRepBuilderAPI_MakeWire that makes sure segments are connected either by moving end points or by adding intermediate segments + class wire_builder { + private: + BRepBuilderAPI_MakeWire mw_; + double p_; + bool override_next_; + gp_Pnt next_override_; + const IfcUtil::IfcBaseClass* inst_; + + public: + wire_builder(double p, const IfcUtil::IfcBaseClass* inst = 0) : p_(p), override_next_(false), inst_(inst) {} + + void operator()(const TopoDS_Shape& a); + + void operator()(const TopoDS_Shape& a, const TopoDS_Shape& b, bool last); + + const TopoDS_Wire& wire() { return mw_.Wire(); } + }; + + template + void shape_pair_enumerate(TopTools_ListIteratorOfListOfShape& it, Fn& fn, bool closed) { + bool is_first = true; + TopoDS_Shape first, previous, current; + for (; it.More(); it.Next(), is_first = false) { + current = it.Value(); + if (is_first) { + first = current; + } else { + fn(previous, current, false); + } + previous = current; + } + if (closed) { + fn(current, first, true); + } else { + fn(current); + } + } + + /* + Below is code to deduce the formula below in SageMath + + | R, b = var('R b') + | + | Bxy = R * cos(b), R * sin(b) + | Cxy = R * cos(b/2), R * sin(b/2) + | + | def dot(v, w): + | return v[0] * w[0] + v[1] * w[1] + | + | def norm(v): + | l = sqrt(v[0]^2 + v[1]^2) + | return v[0] / l, v[1] / l + | + | (R - R*dot(norm(Cxy), norm(Bxy))).full_simplify() + */ + + double deflection_for_approximating_circle(double radius, double param); + + bool create_edge_over_curve_with_log_messages(const Handle_Geom_Curve& crv, const double eps, const gp_Pnt& p1, const gp_Pnt& p2, TopoDS_Edge& result); + } +} + +#endif