From 0b3e9e7c5acacac89c1783cd9b363888e303739c Mon Sep 17 00:00:00 2001 From: Dion Moult Date: Wed, 7 Feb 2024 17:49:12 +1100 Subject: [PATCH] Reuse geom iterator triangulation data, don't use optimal OBBs for speed, verts now use gp_Pnt instead of BVH_Vec3d, and optimise add_triangulated. --- .../IfcGeomRepresentation.cpp | 3 +- src/ifcgeom_schema_agnostic/IfcGeomTree.h | 211 +++++++----------- 2 files changed, 78 insertions(+), 136 deletions(-) diff --git a/src/ifcgeom_schema_agnostic/IfcGeomRepresentation.cpp b/src/ifcgeom_schema_agnostic/IfcGeomRepresentation.cpp index 7672eb381f..343da1a082 100644 --- a/src/ifcgeom_schema_agnostic/IfcGeomRepresentation.cpp +++ b/src/ifcgeom_schema_agnostic/IfcGeomRepresentation.cpp @@ -481,7 +481,8 @@ IfcGeom::Representation::Triangulation::Triangulation(const BRep& shape_model) } } - BRepTools::Clean(s); + // Temporarily commented out so we don't need to triangulate it again for clash detection. + // BRepTools::Clean(s); } } diff --git a/src/ifcgeom_schema_agnostic/IfcGeomTree.h b/src/ifcgeom_schema_agnostic/IfcGeomTree.h index bed524e616..fda31dba82 100644 --- a/src/ifcgeom_schema_agnostic/IfcGeomTree.h +++ b/src/ifcgeom_schema_agnostic/IfcGeomTree.h @@ -47,7 +47,6 @@ #include #include #include -#include #include #include #include @@ -262,7 +261,7 @@ namespace IfcGeom { bool is_point_in_shape( const gp_Pnt& v, const opencascade::handle>& bvh, - const std::unordered_map>& verts, + const std::unordered_map>& verts, // In the case of "touching" rays, let's check again! bool should_check_again = false ) const { @@ -340,14 +339,11 @@ namespace IfcGeom { //std::cout << "Ray hits leaf" << std::endl; // Do ray triangle check. for (int j=bvh->BegPrimitive(i); j<=bvh->EndPrimitive(i); ++j) { - const std::array& v123 = verts.at(j); - const BVH_Vec3d& v1 = v123[0]; - const BVH_Vec3d& v2 = v123[1]; - const BVH_Vec3d& v3 = v123[2]; + const std::array& v123 = verts.at(j); - gp_Vec ta(v1[0], v1[1], v1[2]); - gp_Vec tb(v2[0], v2[1], v2[2]); - gp_Vec tc(v3[0], v3[1], v3[2]); + gp_Vec ta(v123[0].X(), v123[0].Y(), v123[0].Z()); + gp_Vec tb(v123[1].X(), v123[1].Y(), v123[2].Z()); + gp_Vec tc(v123[2].X(), v123[2].Y(), v123[3].Z()); /* std::cout << "ray origin " << ray_origin.X() << " " << ray_origin.Y() << " " << ray_origin.Z() << std::endl; @@ -381,7 +377,7 @@ namespace IfcGeom { const gp_Vec& e2, const opencascade::handle>& bvh, const std::unordered_map& valid_tris, - const std::unordered_map>& verts, + const std::unordered_map>& verts, const std::unordered_map& normals ) const { const gp_Vec& ray_origin = e1; @@ -436,19 +432,16 @@ namespace IfcGeom { if ( ! valid_tris.at(j)) { continue; } - const std::array& v123 = verts.at(j); - const BVH_Vec3d& v1 = v123[0]; - const BVH_Vec3d& v2 = v123[1]; - const BVH_Vec3d& v3 = v123[2]; + const std::array& v123 = verts.at(j); const gp_Vec& normal = normals.at(j); if (std::abs(normal.Dot(ray_vector)) < 1e-3) { continue; // This ray is coplanar to the triangle } - gp_Vec ta(v1[0], v1[1], v1[2]); - gp_Vec tb(v2[0], v2[1], v2[2]); - gp_Vec tc(v3[0], v3[1], v3[2]); + gp_Vec ta(v123[0].X(), v123[0].Y(), v123[0].Z()); + gp_Vec tb(v123[1].X(), v123[1].Y(), v123[1].Z()); + gp_Vec tc(v123[2].X(), v123[2].Y(), v123[2].Z()); double at, au, av; // Do box check first? @@ -835,8 +828,8 @@ namespace IfcGeom { const std::unordered_map& valid_tris_a = valid_tris_.find(tA)->second; const std::unordered_map& valid_tris_b = valid_tris_.find(tB)->second; - const std::unordered_map>& verts_a = verts_.find(tA)->second; - const std::unordered_map>& verts_b = verts_.find(tB)->second; + const std::unordered_map>& verts_a = verts_.find(tA)->second; + const std::unordered_map>& verts_b = verts_.find(tB)->second; const std::unordered_map& normals_a = normals_.find(tA)->second; const std::unordered_map& normals_b = normals_.find(tB)->second; @@ -861,16 +854,12 @@ namespace IfcGeom { continue; } - const std::array& verts = verts_a.at(i); - const BVH_Vec3d& v1 = verts[0]; - const BVH_Vec3d& v2 = verts[1]; - const BVH_Vec3d& v3 = verts[2]; + const std::array& verts = verts_a.at(i); + const gp_Pnt& v1_a_pnt = verts[0]; + const gp_Pnt& v2_a_pnt = verts[1]; + const gp_Pnt& v3_a_pnt = verts[2]; const gp_Vec& normal_a = normals_a.at(i); - const gp_Pnt v1_a_pnt(v1[0], v1[1], v1[2]); - const gp_Pnt v2_a_pnt(v2[0], v2[1], v2[2]); - const gp_Pnt v3_a_pnt(v3[0], v3[1], v3[2]); - const std::array points_a = {v1_a_pnt, v2_a_pnt, v3_a_pnt}; std::vector points_in_b; @@ -949,16 +938,9 @@ namespace IfcGeom { continue; } - const std::array& verts = verts_b.at(j); - const BVH_Vec3d& v1_b = verts[0]; - const BVH_Vec3d& v2_b = verts[1]; - const BVH_Vec3d& v3_b = verts[2]; + const std::array& verts = verts_b.at(j); const gp_Vec& normal_b = normals_b.at(j); - const gp_Pnt v1_b_pnt(v1_b[0], v1_b[1], v1_b[2]); - const gp_Pnt v2_b_pnt(v2_b[0], v2_b[1], v2_b[2]); - const gp_Pnt v3_b_pnt(v3_b[0], v3_b[1], v3_b[2]); - tri_count_++; // We're penetrating _into_ a shape, so don't @@ -968,13 +950,12 @@ namespace IfcGeom { continue; } + gp_Vec ta(verts[0].X(), verts[0].Y(), verts[0].Z()); + gp_Vec tb(verts[1].X(), verts[1].Y(), verts[1].Z()); + gp_Vec tc(verts[2].X(), verts[2].Y(), verts[2].Z()); + for (const auto& v : points_in_b) { gp_Vec ray_origin(v.X(), v.Y(), v.Z()); - gp_Vec point_on_b; - - gp_Vec ta(v1_b[0], v1_b[1], v1_b[2]); - gp_Vec tb(v2_b[0], v2_b[1], v2_b[2]); - gp_Vec tc(v3_b[0], v3_b[1], v3_b[2]); /* std::cout << "POINT IN B " << v.X() << " " << v.Y() << " " << v.Z() << std::endl; @@ -989,13 +970,6 @@ namespace IfcGeom { if (intersectRayTriangle(ray_origin, normal_b, ta, tb, tc, at, au, av, false)) { double current_v_protrusion = at; - /* - // What happens now? - if (current_v_protrusion > max_protrusion) { - continue; - } - */ - // std::cout << "We got a current protrusion " << current_v_protrusion << std::endl; if (current_v_protrusion < v_protrusion) { double aw = 1.0f - au - av; // Barycentric coordinate for ta @@ -1005,7 +979,8 @@ namespace IfcGeom { v_protrusion_point = {v.X(), v.Y(), v.Z()}; v_surface_point = {point_on_b.X(), point_on_b.Y(), point_on_b.Z()}; - if ( ! check_all && v_protrusion > tolerance) { + if (( ! check_all && v_protrusion > tolerance) + || v_protrusion > (max_protrusion - 1e-3)) { clash_types_.push_back(0); protrusion_distances_.push_back(v_protrusion); protrusion_points_.push_back(v_protrusion_point); @@ -1020,7 +995,7 @@ namespace IfcGeom { if (v_protrusion != std::numeric_limits::infinity()) { if (v_protrusion > protrusion) { - std::cout << "New actual protrusion winner of " << v_protrusion << std::endl; + // std::cout << "New actual protrusion winner of " << v_protrusion << std::endl; protrusion = v_protrusion; protrusion_point = v_protrusion_point; surface_point = v_surface_point; @@ -1117,8 +1092,8 @@ namespace IfcGeom { const std::unordered_map& valid_tris_a = valid_tris_.find(tA)->second; const std::unordered_map& valid_tris_b = valid_tris_.find(tB)->second; - const std::unordered_map>& verts_a = verts_.find(tA)->second; - const std::unordered_map>& verts_b = verts_.find(tB)->second; + const std::unordered_map>& verts_a = verts_.find(tA)->second; + const std::unordered_map>& verts_b = verts_.find(tB)->second; const std::unordered_map& normals_a = normals_.find(tA)->second; const std::unordered_map& normals_b = normals_.find(tB)->second; @@ -1131,19 +1106,15 @@ namespace IfcGeom { continue; } - const std::array& verts = verts_a.at(i); - const BVH_Vec3d& v1 = verts[0]; - const BVH_Vec3d& v2 = verts[1]; - const BVH_Vec3d& v3 = verts[2]; + const std::array& verts = verts_a.at(i); + const gp_Pnt& v1_a_pnt = verts[0]; + const gp_Pnt& v2_a_pnt = verts[1]; + const gp_Pnt& v3_a_pnt = verts[2]; const gp_Vec& normal_a = normals_a.at(i); - gp_Pnt v1_a_pnt(v1[0], v1[1], v1[2]); - gp_Pnt v2_a_pnt(v2[0], v2[1], v2[2]); - gp_Pnt v3_a_pnt(v3[0], v3[1], v3[2]); - - std::array t1a = {v1[0], v1[1], v1[2]}; - std::array t1b = {v2[0], v2[1], v2[2]}; - std::array t1c = {v3[0], v3[1], v3[2]}; + std::array t1a = {v1_a_pnt.X(), v1_a_pnt.Y(), v1_a_pnt.Z()}; + std::array t1b = {v2_a_pnt.X(), v2_a_pnt.Y(), v2_a_pnt.Z()}; + std::array t1c = {v3_a_pnt.X(), v3_a_pnt.Y(), v3_a_pnt.Z()}; for (const auto& bvh_b_i : bvh_b_is) { for (int j=bvh_b->BegPrimitive(bvh_b_i); j<=bvh_b->EndPrimitive(bvh_b_i); ++j) { @@ -1151,21 +1122,17 @@ namespace IfcGeom { continue; } - const std::array& verts = verts_b.at(j); - const BVH_Vec3d& v1_b = verts[0]; - const BVH_Vec3d& v2_b = verts[1]; - const BVH_Vec3d& v3_b = verts[2]; + const std::array& verts = verts_b.at(j); + const gp_Pnt& v1_b_pnt = verts[0]; + const gp_Pnt& v2_b_pnt = verts[1]; + const gp_Pnt& v3_b_pnt = verts[2]; const gp_Vec& normal_b = normals_b.at(j); tri_count_++; - gp_Pnt v1_b_pnt(v1_b[0], v1_b[1], v1_b[2]); - gp_Pnt v2_b_pnt(v2_b[0], v2_b[1], v2_b[2]); - gp_Pnt v3_b_pnt(v3_b[0], v3_b[1], v3_b[2]); - - std::array t2a = {v1_b[0], v1_b[1], v1_b[2]}; - std::array t2b = {v2_b[0], v2_b[1], v2_b[2]}; - std::array t2c = {v3_b[0], v3_b[1], v3_b[2]}; + std::array t2a = {v1_b_pnt.X(), v1_b_pnt.Y(), v1_b_pnt.Z()}; + std::array t2b = {v2_b_pnt.X(), v2_b_pnt.Y(), v2_b_pnt.Z()}; + std::array t2c = {v3_b_pnt.X(), v3_b_pnt.Y(), v3_b_pnt.Z()}; // Allow a deviation of 0.25 degrees in coplanarity check if (std::abs(normal_a.Dot(normal_b)) >= 0.99999f) { @@ -1330,8 +1297,8 @@ namespace IfcGeom { return false; } - const std::unordered_map>& verts_a = verts_.find(tA)->second; - const std::unordered_map>& verts_b = verts_.find(tB)->second; + const std::unordered_map>& verts_a = verts_.find(tA)->second; + const std::unordered_map>& verts_b = verts_.find(tB)->second; double min_clearance = std::numeric_limits::infinity(); std::array clearance_point1; @@ -1342,29 +1309,29 @@ namespace IfcGeom { const std::vector& bvh_b_is = pair.second; for (int i=bvh_a->BegPrimitive(bvh_a_i); i<=bvh_a->EndPrimitive(bvh_a_i); ++i) { - const std::array& verts = verts_a.at(i); - const BVH_Vec3d& v1 = verts[0]; - const BVH_Vec3d& v2 = verts[1]; - const BVH_Vec3d& v3 = verts[2]; + const std::array& verts = verts_a.at(i); + const gp_Pnt& v1_a_pnt = verts[0]; + const gp_Pnt& v2_a_pnt = verts[1]; + const gp_Pnt& v3_a_pnt = verts[2]; - const gp_Vec v1_a_vec(v1[0], v1[1], v1[2]); - const gp_Vec v2_a_vec(v2[0], v2[1], v2[2]); - const gp_Vec v3_a_vec(v3[0], v3[1], v3[2]); + const gp_Vec v1_a_vec(v1_a_pnt.X(), v1_a_pnt.Y(), v1_a_pnt.Z()); + const gp_Vec v2_a_vec(v2_a_pnt.X(), v2_a_pnt.Y(), v2_a_pnt.Z()); + const gp_Vec v3_a_vec(v3_a_pnt.X(), v3_a_pnt.Y(), v3_a_pnt.Z()); const std::array p = {v1_a_vec, v2_a_vec, v3_a_vec}; for (const auto& bvh_b_i : bvh_b_is) { for (int j=bvh_b->BegPrimitive(bvh_b_i); j<=bvh_b->EndPrimitive(bvh_b_i); ++j) { - const std::array& verts = verts_b.at(j); - const BVH_Vec3d& v1_b = verts[0]; - const BVH_Vec3d& v2_b = verts[1]; - const BVH_Vec3d& v3_b = verts[2]; + const std::array& verts = verts_b.at(j); + const gp_Pnt& v1_b_pnt = verts[0]; + const gp_Pnt& v2_b_pnt = verts[1]; + const gp_Pnt& v3_b_pnt = verts[2]; tri_count_++; - const gp_Vec v1_b_vec(v1_b[0], v1_b[1], v1_b[2]); - const gp_Vec v2_b_vec(v2_b[0], v2_b[1], v2_b[2]); - const gp_Vec v3_b_vec(v3_b[0], v3_b[1], v3_b[2]); + const gp_Vec v1_b_vec(v1_b_pnt.X(), v1_b_pnt.Y(), v1_b_pnt.Z()); + const gp_Vec v2_b_vec(v2_b_pnt.X(), v2_b_pnt.Y(), v2_b_pnt.Z()); + const gp_Vec v3_b_vec(v3_b_pnt.X(), v3_b_pnt.Y(), v3_b_pnt.Z()); const std::array q = {v1_b_vec, v2_b_vec, v3_b_vec}; @@ -1465,7 +1432,7 @@ namespace IfcGeom { // Note that the original add function is also used elsewhere (e.g. boolean_utils.cpp) // We don't want to randomly add triangulated voids in our // tree, so for now this is a separate function. - BRepMesh_IncrementalMesh(s, 1.e-3, false, 0.5); + // BRepMesh_IncrementalMesh(s, 1.e-3, false, 0.5); Bnd_Box b; BRepBndLib::AddClose(s, b); @@ -1473,15 +1440,14 @@ namespace IfcGeom { shapes_[t] = s; Bnd_OBB obb; - BRepBndLib::AddOBB(s, obb, true, true, false); + // If IsOptimal = True it doubles the execution time. + BRepBndLib::AddOBB(s, obb, true, false, false); obbs_[t] = obb; max_protrusions_[t] = std::min(std::min(obb.XHSize(), obb.YHSize()), obb.ZHSize()) * 2; - BVH_BoxSet* boxset = new BVH_BoxSet(); BRepExtrema_ShapeList shape_list; - std::unordered_map faces; std::unordered_map is_reversed; TopExp_Explorer exp_f; @@ -1489,39 +1455,22 @@ namespace IfcGeom { for (exp_f.Init(s, TopAbs_FACE); exp_f.More(); exp_f.Next()) { shape_list.Append(exp_f.Current()); - Bnd_Box aabb; - BRepBndLib::Add(exp_f.Current(), aabb); - double x, y, z, X, Y, Z; - aabb.Get(x, y, z, X, Y, Z); - const BVH_Box::BVH_VecNt min(x, y, z); - const BVH_Box::BVH_VecNt max(X, Y, Z); - BVH_Box bvhBox(min, max); - boxset->Add(i, bvhBox); - TopoDS_Face f = TopoDS::Face(exp_f.Current()); - faces[i] = f; is_reversed[i] = f.Orientation() == TopAbs_REVERSED; i++; } - /* Option 1: Builder? - BVH_Tree* bvh = new BVH_Tree(); - BVH_Box bvhBox2; // What's the point of this? - BVH_LinearBuilder builder; - builder.Build(boxset, bvh, bvhBox2); - */ - - /* Option 2: Box set works, but ends up still comparing over 17 billion tri pairs - const opencascade::handle>& bvh = boxset->BVH(); - */ - - // Option 3: Triangle set - down to 96 million pairs BRepExtrema_TriangleSet triangle_set(shape_list); const opencascade::handle>& bvh = triangle_set.BVH(); std::unordered_map valid_tris; - std::unordered_map> verts; + std::unordered_map> verts; std::unordered_map normals; + + valid_tris.reserve(triangle_set.Size()); + verts.reserve(triangle_set.Size()); + normals.reserve(triangle_set.Size()); + for (int i=0; i Precision::Confusion()) { + normals[i] = cross_product.Normalized(); valid_tris[i] = true; - } catch (...) { + } else { valid_tris[i] = false; } } @@ -1564,13 +1512,10 @@ namespace IfcGeom { } */ - triangle_sets_[t] = triangle_set; - boxsets_[t] = boxset; bvhs_[t] = bvh; - faces_[t] = faces; - verts_[t] = verts; - normals_[t] = normals; - valid_tris_[t] = valid_tris; + verts_[t] = std::move(verts); + normals_[t] = std::move(normals); + valid_tris_[t] = std::move(valid_tris); } std::vector select_box(const T& t, bool completely_within = false, double extend=-1.e-5) const { @@ -1824,13 +1769,9 @@ namespace IfcGeom { map_t shapes_; std::map obbs_; std::map max_protrusions_; - //std::map*> bvhs_; std::map>> bvhs_; - std::map*> boxsets_; - std::map triangle_sets_; - std::unordered_map> faces_; std::unordered_map> valid_tris_; - std::unordered_map>> verts_; + std::unordered_map>> verts_; std::unordered_map> normals_; bool enable_face_styles_ = false;