diff --git a/src/ifcgeom/ConversionResult.h b/src/ifcgeom/ConversionResult.h index 281b1b3ed2..4babe69519 100644 --- a/src/ifcgeom/ConversionResult.h +++ b/src/ifcgeom/ConversionResult.h @@ -17,8 +17,8 @@ * * ********************************************************************************/ -#ifndef IFCSHAPELIST_H -#define IFCSHAPELIST_H +#ifndef CONVERSIONRESULT_H +#define CONVERSIONRESULT_H #include "../ifcgeom/IfcGeomRenderStyles.h" #include "../ifcgeom/ConversionSettings.h" @@ -27,6 +27,44 @@ #include #include +struct EdgeKey { + int v1, v2; + + // These are not part of the hash or equality, + // but retained to easily created a directed + // graph of the original boundary edges. Since + // the boundary edges are exactly those with + // count=1 we don't need to worry about + // conflicting original vertex indices. + int ov1, ov2; + + EdgeKey(int a, int b) + : ov1(a) + , ov2(b) + { + if (a < b) { + v1 = a; + v2 = b; + } else { + v1 = b; + v2 = a; + } + } + + bool operator==(const EdgeKey& other) const { + return v1 == other.v1 && v2 == other.v2; + } +}; + +namespace std { + template <> + struct hash { + std::size_t operator()(const EdgeKey& ek) const { + return std::hash()(ek.v1) ^ std::hash()(ek.v2); + } + }; +} + namespace IfcGeom { namespace Representation { @@ -296,6 +334,80 @@ namespace IfcGeom { namespace util { // @todo this is now moved to occt kernel, do we need something similar in cgal? // bool flatten_shape_list(const IfcGeom::ConversionResults& shapes, TopoDS_Shape& result, bool fuse, double tol); + + // Function to find boundary loops from triangles + template + std::vector> find_boundary_loops(const std::vector& positions, const std::vector>& triangles) { + std::unordered_map edge_count; + + // Count how many triangles each edge belongs to + for (const auto& triangle : triangles) { + int v1, v2, v3; + std::tie(v1, v2, v3) = triangle; + + edge_count[{v1, v2}]++; + edge_count[{v2, v3}]++; + edge_count[{v3, v1}]++; + } + + // Boundary edges have count 1 + std::vector boundary_edges; + for (auto& p : edge_count) { + if (p.second == 1) { + boundary_edges.push_back(p.first); + } + } + + // We retained original directed edges so we build + // a mapping out of these directed edges. + std::unordered_map vertex_successors; + for (const auto& e : boundary_edges) { + vertex_successors[e.ov1] = e.ov2; + } + + std::vector> loops; + while (!vertex_successors.empty()) { + loops.emplace_back(); + auto it = vertex_successors.begin(); + loops.back() = { it->first, it->second }; + vertex_successors.erase(it); + + int current = loops.back().back(); + while (!vertex_successors.empty() && current != loops.back().front()) { + auto next = vertex_successors[current]; + if (loops.back().front() != next) { + loops.back().push_back(next); + } + vertex_successors.erase(current); + current = next; + } + } + + // Sort the loops by smallest x-coord of their constituent positions + // In order to put the outermost loop in front + if (loops.size() > 1) { + std::vector> min_xs; + for (auto& l : loops) { + NT min_x = std::numeric_limits::infinity(); + for (auto& i : l) { + const auto& x = positions[i * 3]; + if (x < min_x) { + min_x = x; + } + } + min_xs.push_back({ min_x, min_xs.size() }); + } + std::sort(min_xs.begin(), min_xs.end()); + decltype(loops) loops_copy; + for (auto& p : min_xs) { + loops_copy.emplace_back(std::move(loops[p.second])); + } + std::swap(loops, loops_copy); + } + + return loops; + } } } + #endif diff --git a/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp b/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp index 7a93bd00a9..1e438d2cd9 100644 --- a/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp +++ b/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp @@ -19,7 +19,64 @@ using ifcopenshell::geometry::NumberEpeck; #define NumberType NumberEpeck #endif -ifcopenshell::geometry::CgalShape::CgalShape(const cgal_shape_t & shape, bool convex) { +typedef CGAL::Polyhedron_3 Polyhedron; +typedef Polyhedron::Facet_const_handle Facet_const_handle; +typedef Polyhedron::Halfedge_around_facet_const_circulator Halfedge_around_facet_circulator; + +namespace { + bool are_facets_coplanar(const Facet_const_handle& f1, const Facet_const_handle& f2) { + // Function to determine if two facets are coplanar + // You can use the normal vectors and the equation of the planes to determine coplanarity + auto normal_1 = CGAL::normal(f1->halfedge()->vertex()->point(), + f1->halfedge()->next()->vertex()->point(), + f1->halfedge()->next()->next()->vertex()->point()); + + auto normal_2 = CGAL::normal(f2->halfedge()->vertex()->point(), + f2->halfedge()->next()->vertex()->point(), + f2->halfedge()->next()->next()->vertex()->point()); + + return CGAL::collinear(CGAL::ORIGIN + decltype(normal_1)(0., 0., 0.), CGAL::ORIGIN + normal_1, CGAL::ORIGIN + normal_2); + } + + void partition_coplanar_components(const Polyhedron& shape, + std::vector>& components) { + std::set visited; + + for (auto& face : shape.facet_handles()) { + if (visited.find(face) != visited.end()) { + continue; + } + + // Create a new component for coplanar facets + std::set component; + std::queue queue; + + queue.push(face); + visited.insert(face); + + while (!queue.empty()) { + Facet_const_handle current = queue.front(); + queue.pop(); + + component.insert(current); + + // Iterate over neighboring facets + Halfedge_around_facet_circulator he = current->facet_begin(); + do { + Facet_const_handle neighbour = he->opposite()->face(); + if (neighbour != nullptr && visited.find(neighbour) == visited.end() && are_facets_coplanar(current, neighbour)) { + queue.push(neighbour); + visited.insert(neighbour); + } + } while (++he != current->facet_begin()); + } + + components.push_back(component); + } + } +} + +ifcopenshell::geometry::CgalShape::CgalShape(const cgal_shape_t& shape, bool convex) { shape_ = shape; convex_tag_ = convex; @@ -141,6 +198,17 @@ void ifcopenshell::geometry::CgalShape::Triangulate(ifcopenshell::geometry::Sett } + // Facet -> planar component map for determining which + // edges are to be registered. + std::vector> components; + partition_coplanar_components(s, components); + std::map facet_to_component; + for (auto it = components.begin(); it != components.end(); ++it) { + for (auto& f : *it) { + facet_to_component[f] = it; + } + } + // std::map vertex_normals; // boost::associative_property_map> vertex_normals_map(vertex_normals); @@ -161,6 +229,8 @@ void ifcopenshell::geometry::CgalShape::Triangulate(ifcopenshell::geometry::Sett typedef std::tuple postion_normal; std::map welds; + std::set> registered_edges; + int num_faces = 0, num_vertices = 0; for (auto &face : faces(s)) { if (!face->is_triangle()) { @@ -169,6 +239,7 @@ void ifcopenshell::geometry::CgalShape::Triangulate(ifcopenshell::geometry::Sett } CGAL::Polyhedron_3::Halfedge_around_facet_const_circulator current_halfedge = face->facet_begin(); int vertexidx[3]; + bool is_face_boundary[3]; int i = 0; do { postion_normal pn = { @@ -203,13 +274,32 @@ void ifcopenshell::geometry::CgalShape::Triangulate(ifcopenshell::geometry::Sett vidx = it->second; } - vertexidx[i++] = (int) vidx; + vertexidx[i] = (int)vidx; + is_face_boundary[i] = facet_to_component[face] != facet_to_component[current_halfedge->opposite()->face()]; + ++i; ++num_vertices; ++current_halfedge; } while (current_halfedge != face->facet_begin()); t->addFace(item_id, surface_style_id, vertexidx[0], vertexidx[1], vertexidx[2]); + for (size_t i = 0; i < 3; ++i) { + if (is_face_boundary[i]) { + // In CGAL, the vertex of a halfedge is the incident vertex, i.e + // the second vertex of the edge, so in order to get corresponding + // vertex and edge indices we need to find vertexids (i-1, i) for + // the boundary registered in i. + auto a = vertexidx[(i + 2) % 3]; + auto b = vertexidx[(i + 3) % 3]; + if (a > b) { + std::swap(a, b); + } + if (registered_edges.find({ a, b }) == registered_edges.end()) { + registered_edges.insert({ a,b }); + t->registerEdge(item_id, a, b); + } + } + } ++num_faces; } diff --git a/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.cpp b/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.cpp index d5e6e2c61a..6781f50c9b 100644 --- a/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.cpp +++ b/src/ifcgeom/kernels/opencascade/OpenCascadeConversionResult.cpp @@ -31,44 +31,6 @@ using IfcGeom::OpaqueCoordinate; using IfcGeom::NumberNativeDouble; using IfcGeom::ConversionResultShape; -struct EdgeKey { - int v1, v2; - - // These are not part of the hash or equality, - // but retained to easily created a directed - // graph of the original boundary edges. Since - // the boundary edges are exactly those with - // count=1 we don't need to worry about - // conflicting original vertex indices. - int ov1, ov2; - - EdgeKey(int a, int b) - : ov1(a) - , ov2(b) - { - if (a < b) { - v1 = a; - v2 = b; - } else { - v1 = b; - v2 = a; - } - } - - bool operator==(const EdgeKey& other) const { - return v1 == other.v1 && v2 == other.v2; - } -}; - -namespace std { - template <> - struct hash { - std::size_t operator()(const EdgeKey& ek) const { - return std::hash()(ek.v1) ^ std::hash()(ek.v2); - } - }; -} - namespace { // We bypass the conversion to gp_GTrsf, because it does not work void taxonomy_transform(const Eigen::Matrix4d* m, gp_XYZ& xyz) { @@ -80,78 +42,6 @@ namespace { xyz.ChangeData()[2] = v2(2); } } - - // Function to find boundary loops from triangles - std::vector> find_boundary_loops(const std::vector& positions, const std::vector>& triangles) { - std::unordered_map edge_count; - - // Count how many triangles each edge belongs to - for (const auto& triangle : triangles) { - int v1, v2, v3; - std::tie(v1, v2, v3) = triangle; - - edge_count[{v1, v2}]++; - edge_count[{v2, v3}]++; - edge_count[{v3, v1}]++; - } - - // Boundary edges have count 1 - std::vector boundary_edges; - for (auto& p : edge_count) { - if (p.second == 1) { - boundary_edges.push_back(p.first); - } - } - - // We retained original directed edges so we build - // a mapping out of these directed edges. - std::unordered_map vertex_successors; - for (const auto& e : boundary_edges) { - vertex_successors[e.ov1] = e.ov2; - } - - std::vector> loops; - while (!vertex_successors.empty()) { - loops.emplace_back(); - auto it = vertex_successors.begin(); - loops.back() = { it->first, it->second }; - vertex_successors.erase(it); - - int current = loops.back().back(); - while (!vertex_successors.empty() && current != loops.back().front()) { - auto next = vertex_successors[current]; - if (loops.back().front() != next) { - loops.back().push_back(next); - } - vertex_successors.erase(current); - current = next; - } - } - - // Sort the loops by smallest x-coord of their constituent positions - // In order to put the outermost loop in front - if (loops.size() > 1) { - std::vector> min_xs; - for (auto& l : loops) { - double min_x = std::numeric_limits::infinity(); - for (auto& i : l) { - const auto& x = positions[i * 3]; - if (x < min_x) { - min_x = x; - } - } - min_xs.push_back({ min_x, min_xs.size() }); - } - std::sort(min_xs.begin(), min_xs.end()); - decltype(loops) loops_copy; - for (auto& p : min_xs) { - loops_copy.emplace_back(std::move(loops[p.second])); - } - std::swap(loops, loops_copy); - } - - return loops; - } } void ifcopenshell::geometry::OpenCascadeShape::Triangulate(ifcopenshell::geometry::Settings settings, const ifcopenshell::geometry::taxonomy::matrix4& place, IfcGeom::Representation::Triangulation* t, int item_id, int surface_style_id) const { @@ -313,7 +203,7 @@ void ifcopenshell::geometry::OpenCascadeShape::Triangulate(ifcopenshell::geometr } if (polyhedral_output_without_holes || polyhedral_output_with_holes) { - auto loops = find_boundary_loops(t->verts(), triangle_indices); + auto loops = IfcGeom::util::find_boundary_loops(t->verts(), triangle_indices); if (polyhedral_output_without_holes) { if (!loops.empty() && !loops[0].empty()) { t->addFace(item_id, surface_style_id, loops[0]);