diff --git a/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp b/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp index 1e438d2cd9..f6cd691bcc 100644 --- a/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp +++ b/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp @@ -24,6 +24,28 @@ typedef Polyhedron::Facet_const_handle Facet_const_handle; typedef Polyhedron::Halfedge_around_facet_const_circulator Halfedge_around_facet_circulator; namespace { + template + CGAL::Direction_3 newell(Facet& face) { + typename Kernel_::FT a(0), b(0), c(0); + CGAL::Polyhedron_3::Halfedge_around_facet_const_circulator current_halfedge = face.facet_begin(); + do { + auto& curr = current_halfedge->vertex()->point(); + auto& next = current_halfedge->next()->vertex()->point(); + a += (curr.y() - next.y()) * (curr.z() + next.z()); + b += (curr.z() - next.z()) * (curr.x() + next.x()); + c += (curr.x() - next.x()) * (curr.y() + next.y()); + } while (++current_halfedge != face.facet_begin()); + return CGAL::Direction_3(a, b, c); + } + + struct Plane_equation { + template + typename Facet::Plane_3 operator()(Facet& face) { + typename Facet::Halfedge_handle h = face.halfedge(); + return typename Facet::Plane_3(h->vertex()->point(), newell(face)); + } + }; + 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 @@ -80,15 +102,32 @@ ifcopenshell::geometry::CgalShape::CgalShape(const cgal_shape_t& shape, bool con shape_ = shape; convex_tag_ = convex; + std::set faces_to_remove; + for (const auto& face : CGAL::faces(*shape_)) { // @todo O^2 alert! Use aabb tree or box intersections - bool has_self_intersection = false; + + auto V = newell(*face).to_vector(); + if (V.squared_length() == 0) { + Logger::Warning("Removed face due to self-intersections"); + faces_to_remove.insert(face); + continue; + } + auto C = face->halfedge()->vertex()->point(); + auto transform_point = [&V, &C](const auto& p) { + auto dv = p - C; + return C + (dv - (dv * V) * V); + }; + for (auto& he1 : CGAL::halfedges_around_face(face->halfedge(), *shape_)) { CGAL::Segment_3 s1; { const auto& source = he1->vertex()->point(); const auto& target = he1->next()->vertex()->point(); - s1 = { source, target }; + s1 = { + transform_point(source), + transform_point(target) + }; } for (auto& he2 : CGAL::halfedges_around_face(face->halfedge(), *shape_)) { if (he1 == he2 || he1->next() == he2 || he2->next() == he1) { @@ -99,17 +138,27 @@ ifcopenshell::geometry::CgalShape::CgalShape(const cgal_shape_t& shape, bool con { const auto& source = he2->vertex()->point(); const auto& target = he2->next()->vertex()->point(); - s2 = { source, target }; + s2 = { + transform_point(source), + transform_point(target) + }; } - if (CGAL::do_intersect(s1, s2)) { - has_self_intersection = true; - break; + if (CGAL::to_double((s1.start() - s2.end()).squared_length()) < 1.e-12 || + CGAL::to_double((s1.end() - s2.end()).squared_length()) < 1.e-12 || + CGAL::to_double((s1.start() - s2.start()).squared_length()) < 1.e-12 || + CGAL::to_double((s1.end() - s2.start()).squared_length()) < 1.e-12 || + CGAL::do_intersect(s1, s2)) + { + Logger::Warning("Removed face due to self-intersections"); + faces_to_remove.insert(face); } } } + } - if (has_self_intersection) { - throw std::runtime_error("Self-intersection in facet boundary, not attempting triangulation"); + { + for (auto& face : faces_to_remove) { + CGAL::Euler::remove_face(face->halfedge(), *shape_); } } @@ -444,30 +493,6 @@ OpaqueCoordinate<3> ifcopenshell::geometry::CgalShape::position() } } -namespace { - template - CGAL::Direction_3 newell(Facet& face) { - typename Kernel_::FT a(0), b(0), c(0); - CGAL::Polyhedron_3::Halfedge_around_facet_const_circulator current_halfedge = face.facet_begin(); - do { - auto& curr = current_halfedge->vertex()->point(); - auto& next = current_halfedge->next()->vertex()->point(); - a += (curr.y() - next.y()) * (curr.z() + next.z()); - b += (curr.z() - next.z()) * (curr.x() + next.x()); - c += (curr.x() - next.x()) * (curr.y() + next.y()); - } while (++current_halfedge != face.facet_begin()); - return CGAL::Direction_3(a, b, c); - } - - struct Plane_equation { - template - typename Facet::Plane_3 operator()(Facet& face) { - typename Facet::Halfedge_handle h = face.halfedge(); - return typename Facet::Plane_3(h->vertex()->point(), newell(face)); - } - }; -} - OpaqueCoordinate<3> ifcopenshell::geometry::CgalShape::axis() { to_poly();