From 2cb280932c61aeab212dca51b5c15006f622284b Mon Sep 17 00:00:00 2001 From: Dion Moult Date: Thu, 8 Feb 2024 14:36:20 +1100 Subject: [PATCH] Indirection using vertex indices saves 20% memory but speed is now worse than when we used unordered maps. --- src/ifcgeom_schema_agnostic/IfcGeomTree.h | 165 ++++++++++++---------- 1 file changed, 90 insertions(+), 75 deletions(-) diff --git a/src/ifcgeom_schema_agnostic/IfcGeomTree.h b/src/ifcgeom_schema_agnostic/IfcGeomTree.h index d32f266488..36333cb8d1 100644 --- a/src/ifcgeom_schema_agnostic/IfcGeomTree.h +++ b/src/ifcgeom_schema_agnostic/IfcGeomTree.h @@ -139,7 +139,8 @@ namespace IfcGeom { bool is_point_in_shape( const gp_Pnt& v, const opencascade::handle>& bvh, - const std::vector>& verts, + const std::vector>& tris, + const std::vector& verts, // In the case of "touching" rays, let's check again! bool should_check_again = false ) const { @@ -217,11 +218,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[j]; + const std::array& tri = tris[j]; - 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()); + gp_Vec ta(verts[tri[0]].XYZ()); + gp_Vec tb(verts[tri[1]].XYZ()); + gp_Vec tc(verts[tri[2]].XYZ()); /* std::cout << "ray origin " << ray_origin.X() << " " << ray_origin.Y() << " " << ray_origin.Z() << std::endl; @@ -255,7 +256,8 @@ namespace IfcGeom { const gp_Vec& e2, const opencascade::handle>& bvh, const std::vector& valid_tris, - const std::vector>& verts, + const std::vector>& tris, + const std::vector& verts, const std::vector& normals ) const { const gp_Vec& ray_origin = e1; @@ -310,16 +312,16 @@ namespace IfcGeom { if ( ! valid_tris[j]) { continue; } - const std::array& v123 = verts[j]; + const std::array& tri = tris[j]; const gp_Vec& normal = normals[j]; if (std::abs(normal.Dot(ray_vector)) < 1e-3) { continue; // This ray is coplanar to the triangle } - 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()); + gp_Vec ta(verts[tri[0]].XYZ()); + gp_Vec tb(verts[tri[1]].XYZ()); + gp_Vec tc(verts[tri[2]].XYZ()); double at, au, av; // Do box check first? @@ -478,14 +480,16 @@ namespace IfcGeom { const std::vector& valid_tris_a = valid_tris_.find(tA)->second; const std::vector& valid_tris_b = valid_tris_.find(tB)->second; - const std::vector>& verts_a = verts_.find(tA)->second; - const std::vector>& verts_b = verts_.find(tB)->second; + const std::vector>& tris_a = tris_.find(tA)->second; + const std::vector>& tris_b = tris_.find(tB)->second; + const std::vector& verts_a = verts_.find(tA)->second; + const std::vector& verts_b = verts_.find(tB)->second; const std::vector& normals_a = normals_.find(tA)->second; const std::vector& normals_b = normals_.find(tB)->second; // ~10% faster? - std::unordered_set points_in_b_cache; - std::unordered_set points_not_in_b_cache; + std::unordered_set points_in_b_cache; + std::unordered_set points_not_in_b_cache; double protrusion = -std::numeric_limits::infinity(); std::array protrusion_point; @@ -504,49 +508,41 @@ namespace IfcGeom { continue; } - const std::array& verts = verts_a[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[i]; - - const std::array points_a = {v1_a_pnt, v2_a_pnt, v3_a_pnt}; + const std::array& tri = tris_a[i]; std::vector points_in_b; - for (const auto& v : points_a) { - if (points_not_in_b_cache.find(v) != points_not_in_b_cache.end()) { + for (int v_id : tri) { + if (points_not_in_b_cache.find(v_id) != points_not_in_b_cache.end()) { continue; } - if (points_in_b_cache.find(v) != points_in_b_cache.end()) { + const gp_Pnt& v = verts_a[v_id]; + + if (points_in_b_cache.find(v_id) != points_in_b_cache.end()) { points_in_b.push_back(v); continue; } if (obb_b.IsOut(v)) { - points_not_in_b_cache.insert(v); + points_not_in_b_cache.insert(v_id); continue; } - if (is_point_in_shape(v, bvh_b, verts_b) - && is_point_in_shape(v, bvh_b, verts_b, true)) { + if (is_point_in_shape(v, bvh_b, tris_b, verts_b) + && is_point_in_shape(v, bvh_b, tris_b, verts_b, true)) { points_in_b.push_back(v); - points_in_b_cache.insert(v); + points_in_b_cache.insert(v_id); } else { - points_not_in_b_cache.insert(v); + points_not_in_b_cache.insert(v_id); } } - gp_Vec v1_a_vec(v1_a_pnt.X(), v1_a_pnt.Y(), v1_a_pnt.Z()); - gp_Vec v2_a_vec(v2_a_pnt.X(), v2_a_pnt.Y(), v2_a_pnt.Z()); - gp_Vec v3_a_vec(v3_a_pnt.X(), v3_a_pnt.Y(), v3_a_pnt.Z()); - - double v_protrusion = std::numeric_limits::infinity(); - std::array v_protrusion_point; - std::array v_surface_point; - // If there are no points in b, this may be a "piercing" triangle. if (points_in_b.empty()) { + gp_Vec v1_a_vec(verts_a[tri[0]].XYZ()); + gp_Vec v2_a_vec(verts_a[tri[1]].XYZ()); + gp_Vec v3_a_vec(verts_a[tri[2]].XYZ()); + // Protrusions take priority over piercings. We only check for piercings if: // - This is a piercing triangle (e.g. no points in b) // - No protrusion was already found @@ -555,9 +551,9 @@ namespace IfcGeom { std::array< std::tuple, std::array>, 3 > pierce_results = { - pierce_shape(v1_a_vec, v2_a_vec, bvh_b, valid_tris_b, verts_b, normals_b), - pierce_shape(v1_a_vec, v3_a_vec, bvh_b, valid_tris_b, verts_b, normals_b), - pierce_shape(v2_a_vec, v3_a_vec, bvh_b, valid_tris_b, verts_b, normals_b) + pierce_shape(v1_a_vec, v2_a_vec, bvh_b, valid_tris_b, tris_b, verts_b, normals_b), + pierce_shape(v1_a_vec, v3_a_vec, bvh_b, valid_tris_b, tris_b, verts_b, normals_b), + pierce_shape(v2_a_vec, v3_a_vec, bvh_b, valid_tris_b, tris_b, verts_b, normals_b) }; for (const auto& [p_dist, p_min, p_max] : pierce_results) { @@ -581,6 +577,11 @@ namespace IfcGeom { continue; } + const gp_Vec& normal_a = normals_a[i]; + double v_protrusion = std::numeric_limits::infinity(); + std::array v_protrusion_point; + std::array v_surface_point; + // Check for protrusions. 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) { @@ -588,7 +589,7 @@ namespace IfcGeom { continue; } - const std::array& verts = verts_b[j]; + const std::array& tri = tris_b[j]; const gp_Vec& normal_b = normals_b[j]; tri_count_++; @@ -600,12 +601,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()); + gp_Vec ta(verts_b[tri[0]].XYZ()); + gp_Vec tb(verts_b[tri[1]].XYZ()); + gp_Vec tc(verts_b[tri[2]].XYZ()); for (const auto& v : points_in_b) { - gp_Vec ray_origin(v.X(), v.Y(), v.Z()); + gp_Vec ray_origin(v.XYZ()); /* std::cout << "POINT IN B " << v.X() << " " << v.Y() << " " << v.Z() << std::endl; @@ -694,8 +695,10 @@ namespace IfcGeom { const std::vector& valid_tris_a = valid_tris_.find(tA)->second; const std::vector& valid_tris_b = valid_tris_.find(tB)->second; - const std::vector>& verts_a = verts_.find(tA)->second; - const std::vector>& verts_b = verts_.find(tB)->second; + const std::vector>& tris_a = tris_.find(tA)->second; + const std::vector>& tris_b = tris_.find(tB)->second; + const std::vector& verts_a = verts_.find(tA)->second; + const std::vector& verts_b = verts_.find(tB)->second; const std::vector& normals_a = normals_.find(tA)->second; const std::vector& normals_b = normals_.find(tB)->second; @@ -708,10 +711,10 @@ namespace IfcGeom { continue; } - const std::array& verts = verts_a[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 std::array& tri = tris_a[i]; + const gp_Pnt& v1_a_pnt = verts_a[tri[0]]; + const gp_Pnt& v2_a_pnt = verts_a[tri[1]]; + const gp_Pnt& v3_a_pnt = verts_a[tri[2]]; const gp_Vec& normal_a = normals_a[i]; const gp_Vec v1_a_vec(v1_a_pnt.XYZ()); @@ -724,10 +727,10 @@ namespace IfcGeom { continue; } - const std::array& verts = verts_b[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 std::array& tri = tris_b[j]; + const gp_Pnt& v1_b_pnt = verts_b[tri[0]]; + const gp_Pnt& v2_b_pnt = verts_b[tri[1]]; + const gp_Pnt& v3_b_pnt = verts_b[tri[2]]; const gp_Vec& normal_b = normals_b[j]; tri_count_++; @@ -838,8 +841,10 @@ namespace IfcGeom { return false; } - const std::vector>& verts_a = verts_.find(tA)->second; - const std::vector>& verts_b = verts_.find(tB)->second; + const std::vector>& tris_a = tris_.find(tA)->second; + const std::vector>& tris_b = tris_.find(tB)->second; + const std::vector& verts_a = verts_.find(tA)->second; + const std::vector& verts_b = verts_.find(tB)->second; double min_clearance = std::numeric_limits::infinity(); std::array clearance_point1; @@ -850,10 +855,10 @@ 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[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 std::array& tri = tris_a[i]; + const gp_Pnt& v1_a_pnt = verts_a[tri[0]]; + const gp_Pnt& v2_a_pnt = verts_a[tri[1]]; + const gp_Pnt& v3_a_pnt = verts_a[tri[2]]; const gp_Vec v1_a_vec(v1_a_pnt.XYZ()); const gp_Vec v2_a_vec(v2_a_pnt.XYZ()); @@ -863,10 +868,10 @@ namespace IfcGeom { 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[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 std::array& tri = tris_b[j]; + const gp_Pnt& v1_b_pnt = verts_b[tri[0]]; + const gp_Pnt& v2_b_pnt = verts_b[tri[1]]; + const gp_Pnt& v3_b_pnt = verts_b[tri[2]]; tri_count_++; @@ -1003,25 +1008,33 @@ namespace IfcGeom { const opencascade::handle>& bvh = triangle_set.BVH(); std::vector valid_tris(triangle_set.Size(), true); - std::vector> verts(triangle_set.Size()); + std::vector> tris(triangle_set.Size()); std::vector normals(triangle_set.Size()); + const BVH_Array3d& vertices = triangle_set.GetVertices(); + std::vector verts(vertices.size()); + int i = 0; + for (const auto& v : vertices) { + verts[i] = gp_Pnt(v[0], v[1], v[2]); + i++; + } + for (int i=0; i indices; + triangle_set.GetVtxIndices(i, indices); if (is_reversed[triangle_set.GetFaceID(i)]) { - triangle_set.GetVertices(i, v1, v3, v2); - } else { - triangle_set.GetVertices(i, v1, v2, v3); + std::swap(indices[0], indices[2]); } - gp_Pnt v1_pnt(v1[0], v1[1], v1[2]); - gp_Pnt v2_pnt(v2[0], v2[1], v2[2]); - gp_Pnt v3_pnt(v3[0], v3[1], v3[2]); - verts[i] = {v1_pnt, v2_pnt, v3_pnt}; + + gp_Pnt& v1_pnt = verts[indices[0]]; + gp_Pnt& v2_pnt = verts[indices[1]]; + gp_Pnt& v3_pnt = verts[indices[2]]; gp_Vec dir1(v1_pnt, v2_pnt); gp_Vec dir2(v1_pnt, v3_pnt); gp_Vec cross_product = dir1.Crossed(dir2); if (cross_product.Magnitude() > Precision::Confusion()) { normals[i] = cross_product.Normalized(); + tris[i] = {indices[0], indices[1], indices[2]}; } else { valid_tris[i] = false; } @@ -1047,6 +1060,7 @@ namespace IfcGeom { */ bvhs_[t] = bvh; + tris_[t] = std::move(tris); verts_[t] = std::move(verts); normals_[t] = std::move(normals); valid_tris_[t] = std::move(valid_tris); @@ -1305,7 +1319,8 @@ namespace IfcGeom { std::map max_protrusions_; std::map>> bvhs_; std::unordered_map> valid_tris_; - std::unordered_map>> verts_; + std::unordered_map>> tris_; + std::unordered_map> verts_; std::unordered_map> normals_; bool enable_face_styles_ = false;