From 20865d2544041b23af9ced5f41968873659bcd77 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Fri, 26 Feb 2021 15:40:46 +0100 Subject: [PATCH 1/4] robustify side_of test --- .../Corefinement/Face_graph_output_builder.h | 61 +++++++++++------ .../internal/Corefinement/face_graph_utils.h | 68 +++++++++++++++++++ 2 files changed, 108 insertions(+), 21 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h index 6f57675de59..306f0d9144a 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h @@ -1109,7 +1109,14 @@ public: CGAL::Bounded_side in_tm2 = is_tm2_inside_out ? ON_UNBOUNDED_SIDE : ON_BOUNDED_SIDE; - Side_of_triangle_mesh inside_tm2(tm2, vpm2); + typedef typename Nodes_vector::Exact_kernel Exact_kernel; + typedef Side_of_vpm_helper VPM_helper; + typedef typename VPM_helper::type SOTM_vpm2; + + SOTM_vpm2 sotm_vpm2 = VPM_helper::get_vpm(vertex_to_node_id2, vpm2, nodes); + Side_of_triangle_mesh inside_tm2(tm2, sotm_vpm2); for(face_descriptor f : faces(tm1)) { @@ -1120,31 +1127,36 @@ public: patch_status_not_set_tm1.reset( patch_id ); halfedge_descriptor h = halfedge(f, tm1); Node_id index_p1 = get_node_id(target(h, tm1), vertex_to_node_id1); + std::array fnids = { index_p1, index_p1, index_p1 }; if (index_p1 != NID) { h=next(h, tm1); index_p1 = get_node_id(target(h, tm1), vertex_to_node_id1); + fnids[1]=index_p1; if (index_p1 != NID) { h=next(h, tm1); index_p1 = get_node_id(target(h, tm1), vertex_to_node_id1); + fnids[2]=index_p1; } } + if (index_p1 != NID) { - if (tm1_coplanar_faces.test(f_id)) + if (tm1_coplanar_faces.test(f_id)) // TODO delay the building of the tree? { coplanar_patches_of_tm1.set(patch_id); coplanar_patches_of_tm1_for_union_and_intersection.set(patch_id); } else { - // triangle which is tangent at its 3 vertices - // \todo improve this part which is not robust with a kernel - // with inexact constructions. - Bounded_side position = inside_tm2(centroid(get(vpm1, source(h, tm1)), - get(vpm1, target(h, tm1)), - get(vpm1, target(next(h, tm1), tm1)) )); + typename Exact_kernel::Point_3 e_centroid = + centroid(nodes.exact_node(fnids[0]), + nodes.exact_node(fnids[1]), + nodes.exact_node(fnids[2])); + + Bounded_side position = inside_tm2(e_centroid); + CGAL_assertion( position != ON_BOUNDARY); if ( position == in_tm2 ) is_patch_inside_tm2.set(patch_id); @@ -1152,9 +1164,7 @@ public: } else { - // TODO: tm2 might have been modified and an inexact vpm will - // provide a non-robust result. - Bounded_side position = inside_tm2( get(vpm1, target(h, tm1))); + Bounded_side position = inside_tm2( nodes.to_exact(get(vpm1, target(h, tm1)))); CGAL_assertion( position != ON_BOUNDARY); if ( position == in_tm2 ) is_patch_inside_tm2.set(patch_id); @@ -1171,7 +1181,15 @@ public: CGAL::Bounded_side in_tm1 = is_tm1_inside_out ? ON_UNBOUNDED_SIDE : ON_BOUNDED_SIDE; - Side_of_triangle_mesh inside_tm1(tm1, vpm1); + typedef typename Nodes_vector::Exact_kernel Exact_kernel; + typedef Side_of_vpm_helper VPM_helper; + typedef typename VPM_helper::type SOTM_vpm1; + + SOTM_vpm1 sotm_vpm1 = VPM_helper::get_vpm(vertex_to_node_id1, vpm1, nodes); + Side_of_triangle_mesh inside_tm1(tm1, sotm_vpm1); + for(face_descriptor f : faces(tm2)) { const std::size_t f_id = get(fids2, f); @@ -1181,30 +1199,33 @@ public: patch_status_not_set_tm2.reset( patch_id ); halfedge_descriptor h = halfedge(f, tm2); Node_id index_p2 = get_node_id(target(h, tm2), vertex_to_node_id2); + std::array fnids = { index_p2, index_p2, index_p2 }; if (index_p2 != NID) { h=next(h, tm2); index_p2 = get_node_id(target(h, tm2), vertex_to_node_id2); + fnids[1]=index_p2; if (index_p2 != NID) { h=next(h, tm2); index_p2 = get_node_id(target(h, tm2), vertex_to_node_id2); + fnids[2]=index_p2; } } if (index_p2 != NID) { - if (tm2_coplanar_faces.test(f_id)) + if (tm2_coplanar_faces.test(f_id)) // TODO delay the building of the tree? { coplanar_patches_of_tm2.set(patch_id); coplanar_patches_of_tm2_for_union_and_intersection.set(patch_id); } else { - // triangle which is tangent at its 3 vertices - // \todo improve this part which is not robust with a kernel - // with inexact constructions. - Bounded_side position = inside_tm1(midpoint(get(vpm2, source(h, tm2)), - get(vpm2, target(h, tm2)) )); + typename Exact_kernel::Point_3 e_centroid = + centroid(nodes.exact_node(fnids[0]), + nodes.exact_node(fnids[1]), + nodes.exact_node(fnids[2])); + Bounded_side position = inside_tm1(e_centroid); CGAL_assertion( position != ON_BOUNDARY); if ( position == in_tm1 ) is_patch_inside_tm1.set(patch_id); @@ -1212,9 +1233,7 @@ public: } else { - // TODO: tm1 might have been modified and an inexact vpm will - // provide a non-robust result. - Bounded_side position = inside_tm1( get(vpm2, target(h, tm2))); + Bounded_side position = inside_tm1( nodes.to_exact(get(vpm2, target(h, tm2)))); CGAL_assertion( position != ON_BOUNDARY); if ( position == in_tm1 ) is_patch_inside_tm1.set(patch_id); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h index 4103dcd41ae..c4c07117619 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h @@ -146,6 +146,74 @@ void copy_edge_mark(G&, No_mark&) {} // nothing to do +template +struct Node_vector_exact_vertex_point_map +{ +// map type definitions + typedef typename boost::property_traits::key_type key_type; + typedef typename NodeVector::Exact_kernel Exact_kernel; + typedef typename Exact_kernel::Point_3 value_type; + typedef value_type reference; + typedef boost::readable_property_map_tag category; +// internal type definitions + typedef std::size_t Node_id; + + Node_vector_exact_vertex_point_map(){} + + Node_vector_exact_vertex_point_map(const Node_id_map& node_ids, + const VertexPointMap& vpm, + const NodeVector& node_vector) + : node_ids(&node_ids) + , vpm(&vpm) + , node_vector(&node_vector) + {} + + friend value_type get(Node_vector_exact_vertex_point_map m, key_type k) + { + typename Node_id_map::const_iterator it = m.node_ids->find(k); + if (it == m.node_ids->end()) + return m.node_vector->to_exact( get(*(m.vpm), k) ); + return m.node_vector->exact_node(it->second); + } + + const Node_id_map* node_ids; + const VertexPointMap* vpm; + const NodeVector* node_vector; +}; + +template +struct Side_of_vpm_helper +{ + typedef Node_vector_exact_vertex_point_map type; + static + type get_vpm(const Node_id_map& node_ids, + const VertexPointMap& vpm, + const NodeVector& node_vector) + { + return type(node_ids, vpm, node_vector); + } +}; + +template +struct Side_of_vpm_helper +{ + typedef VertexPointMap type; + static + type get_vpm(const Node_id_map&, + const VertexPointMap& vpm, + const NodeVector&) + { + return vpm; + } +}; + // Parts to get default property maps for output meshes based on the value type // of input vertex point maps. template From 313f682b8c01e8e1648537c550db28ba40b3e806 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 2 Mar 2021 16:16:07 +0100 Subject: [PATCH 2/4] use custom functor for bbox and split primitives aabb tree build runtime is similar to using EPICK --- AABB_tree/include/CGAL/AABB_traits.h | 12 +- .../Corefinement/Face_graph_output_builder.h | 29 ++-- .../internal/Corefinement/face_graph_utils.h | 155 ++++++++++++++++-- 3 files changed, 167 insertions(+), 29 deletions(-) diff --git a/AABB_tree/include/CGAL/AABB_traits.h b/AABB_tree/include/CGAL/AABB_traits.h index 613d08c174b..e443e08a522 100644 --- a/AABB_tree/include/CGAL/AABB_traits.h +++ b/AABB_tree/include/CGAL/AABB_traits.h @@ -428,6 +428,11 @@ public: Closest_point closest_point_object() const {return Closest_point(*this);} Compare_distance compare_distance_object() const {return Compare_distance();} + typedef enum { CGAL_AXIS_X = 0, + CGAL_AXIS_Y = 1, + CGAL_AXIS_Z = 2} Axis; + + static Axis longest_axis(const Bounding_box& bbox); private: /** @@ -446,13 +451,6 @@ private: return internal::Primitive_helper::get_datum(pr,*this).bbox(); } - - typedef enum { CGAL_AXIS_X = 0, - CGAL_AXIS_Y = 1, - CGAL_AXIS_Z = 2} Axis; - - static Axis longest_axis(const Bounding_box& bbox); - /// Comparison functions static bool less_x(const Primitive& pr1, const Primitive& pr2,const AABB_traits& traits) { diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h index 306f0d9144a..2dd5efc3790 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h @@ -1110,13 +1110,17 @@ public: ? ON_UNBOUNDED_SIDE : ON_BOUNDED_SIDE; typedef typename Nodes_vector::Exact_kernel Exact_kernel; - typedef Side_of_vpm_helper VPM_helper; - typedef typename VPM_helper::type SOTM_vpm2; + typedef Side_of_helper VPM_helper; + typedef typename VPM_helper::VPM SOTM_vpm2; + typedef typename VPM_helper::Tree_type Tree_type; SOTM_vpm2 sotm_vpm2 = VPM_helper::get_vpm(vertex_to_node_id2, vpm2, nodes); - Side_of_triangle_mesh inside_tm2(tm2, sotm_vpm2); + Tree_type tree; + VPM_helper::build_tree(tm2, tree, vertex_to_node_id2, fids2, vpm2, nodes); + Side_of_triangle_mesh inside_tm2(tree); for(face_descriptor f : faces(tm1)) { @@ -1182,13 +1186,16 @@ public: ? ON_UNBOUNDED_SIDE : ON_BOUNDED_SIDE; typedef typename Nodes_vector::Exact_kernel Exact_kernel; - typedef Side_of_vpm_helper VPM_helper; - typedef typename VPM_helper::type SOTM_vpm1; + typedef Side_of_helper VPM_helper; + typedef typename VPM_helper::VPM SOTM_vpm1; + typedef typename VPM_helper::Tree_type Tree_type; - SOTM_vpm1 sotm_vpm1 = VPM_helper::get_vpm(vertex_to_node_id1, vpm1, nodes); - Side_of_triangle_mesh inside_tm1(tm1, sotm_vpm1); + Tree_type tree; + VPM_helper::build_tree(tm1, tree, vertex_to_node_id1, fids1, vpm1, nodes); + Side_of_triangle_mesh inside_tm1(tree); for(face_descriptor f : faces(tm2)) { diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h index c4c07117619..fcbd22db729 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h @@ -146,6 +146,8 @@ void copy_edge_mark(G&, No_mark&) {} // nothing to do + +// For exact side_of_triangle_mesh template @@ -183,35 +185,166 @@ struct Node_vector_exact_vertex_point_map const NodeVector* node_vector; }; -template +struct Split_primitives +{ + Split_primitives(TriangleMesh& tm, PPM ppm) + : tm(tm) + , ppm(ppm) + {} + + template + void operator()(PrimitiveIterator first, + PrimitiveIterator beyond, + const CGAL::Bbox_3& bbox) const + { + typedef typename std::iterator_traits::value_type Prmtv; + PrimitiveIterator middle = first + (beyond - first)/2; + typedef typename std::iterator_traits::value_type Prmtv; + switch(TreeTraits::longest_axis(bbox)) + { + case TreeTraits::CGAL_AXIS_X: // sort along x + std::nth_element(first, middle, beyond, [this](const Prmtv& p1, const Prmtv& p2){ return get(ppm, p1.id()).x() < get(ppm,p2.id()).x(); }); + break; + case TreeTraits::CGAL_AXIS_Y: // sort along y + std::nth_element(first, middle, beyond, [this](const Prmtv& p1, const Prmtv& p2){ return get(ppm, p1.id()).y() < get(ppm,p2.id()).y(); }); + break; + case TreeTraits::CGAL_AXIS_Z: // sort along z + std::nth_element(first, middle, beyond, [this](const Prmtv& p1, const Prmtv& p2){ return get(ppm, p1.id()).z() < get(ppm,p2.id()).z(); }); + break; + default: + CGAL_error(); + } + } + TriangleMesh& tm; + PPM ppm; +}; + +// For exact side_of_triangle_mesh +template +struct Compute_bbox { + Compute_bbox(const BPM& bpm) + : bpm(bpm) + {} + + template + CGAL::Bbox_3 operator()(ConstPrimitiveIterator first, + ConstPrimitiveIterator beyond) const + { + CGAL::Bbox_3 bbox = get(bpm, first->id()); + for(++first; first != beyond; ++first) + { + bbox += get(bpm, first->id()); + } + return bbox; + } + BPM bpm; +}; + +// For exact side_of_triangle_mesh +template -struct Side_of_vpm_helper +struct Side_of_helper { - typedef Node_vector_exact_vertex_point_map type; + typedef Node_vector_exact_vertex_point_map VPM; + + typedef CGAL::AABB_face_graph_triangle_primitive Primitive; + typedef CGAL::AABB_traits Traits; + typedef CGAL::AABB_tree Tree_type; + static - type get_vpm(const Node_id_map& node_ids, + VPM get_vpm(const Node_id_map& node_ids, const VertexPointMap& vpm, const NodeVector& node_vector) { - return type(node_ids, vpm, node_vector); + return VPM(node_ids, vpm, node_vector); + } + + template + static + void build_tree(TriangleMesh& tm, + Tree_type& tree, + const Node_id_map& node_ids, + FaceIdMap fid, + const VertexPointMap& vpm, + const NodeVector& node_vector) + { + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + // add primitives + tree.insert(faces(tm).begin(), faces(tm).end(), tm, get_vpm(node_ids, vpm, node_vector)); + + // pre-build bboxes (using approximation) + std::vector face_bboxes(num_faces(tm)); + + auto get_v_box = [&node_ids, &node_vector, &vpm](vertex_descriptor v) + { + typename Node_id_map::const_iterator it = node_ids.find(v); + if (it == node_ids.end()) + return get(vpm, v).bbox(); + return approx(node_vector.exact_node(it->second)).bbox(); + }; + + for (face_descriptor f : faces(tm)) + { + halfedge_descriptor h = halfedge(f, tm); + face_bboxes[get(fid, f)] = get_v_box( source(h, tm) ) + + get_v_box( target(h, tm) ) + + get_v_box( target(next(h, tm), tm) ); + } + + typedef CGAL::Pointer_property_map::type Id_to_box; + Id_to_box id_to_box = CGAL::make_property_map(face_bboxes); + typedef Property_map_binder BPM; + BPM bpm(fid, id_to_box); + Compute_bbox compute_bbox(bpm); + + typedef One_point_from_face_descriptor_map PPM; + PPM ppm(&tm, vpm); + + Split_primitives split_primitives(tm, ppm); + tree.custom_build(compute_bbox, split_primitives); } }; -template -struct Side_of_vpm_helper +struct Side_of_helper { - typedef VertexPointMap type; + typedef VertexPointMap VPM; static - type get_vpm(const Node_id_map&, - const VertexPointMap& vpm, - const NodeVector&) + VPM get_vpm(const Node_id_map&, + const VertexPointMap& vpm, + const NodeVector&) { return vpm; } + + typedef CGAL::AABB_face_graph_triangle_primitive Primitive; + typedef CGAL::AABB_traits Traits; + typedef CGAL::AABB_tree Tree_type; + + + template + static + void build_tree(TriangleMesh& tm, + Tree_type& tree, + const Node_id_map& /* node_ids */, + FaceIdMap /* fid */, + const VertexPointMap& vpm, + const NodeVector& /* node_vector */) + { + tree.insert(faces(tm).begin(), faces(tm).end(), tm, vpm); + tree.build(); + } }; // Parts to get default property maps for output meshes based on the value type From c1dcaac0201de72493579881152e2118c40cc45e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 2 Mar 2021 16:21:06 +0100 Subject: [PATCH 3/4] remove TODOs --- .../internal/Corefinement/Face_graph_output_builder.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h index 2dd5efc3790..341bd670405 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h @@ -1147,7 +1147,7 @@ public: if (index_p1 != NID) { - if (tm1_coplanar_faces.test(f_id)) // TODO delay the building of the tree? + if (tm1_coplanar_faces.test(f_id)) { coplanar_patches_of_tm1.set(patch_id); coplanar_patches_of_tm1_for_union_and_intersection.set(patch_id); @@ -1221,7 +1221,7 @@ public: } if (index_p2 != NID) { - if (tm2_coplanar_faces.test(f_id)) // TODO delay the building of the tree? + if (tm2_coplanar_faces.test(f_id)) { coplanar_patches_of_tm2.set(patch_id); coplanar_patches_of_tm2_for_union_and_intersection.set(patch_id); From fa47c352711a016a299bc990e31d98c16cc2e88e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 3 Mar 2021 08:28:28 +0100 Subject: [PATCH 4/4] remove unused variable --- .../internal/Corefinement/Face_graph_output_builder.h | 1 - 1 file changed, 1 deletion(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h index 341bd670405..3b3d9f82699 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h @@ -1117,7 +1117,6 @@ public: typedef typename VPM_helper::VPM SOTM_vpm2; typedef typename VPM_helper::Tree_type Tree_type; - SOTM_vpm2 sotm_vpm2 = VPM_helper::get_vpm(vertex_to_node_id2, vpm2, nodes); Tree_type tree; VPM_helper::build_tree(tm2, tree, vertex_to_node_id2, fids2, vpm2, nodes); Side_of_triangle_mesh inside_tm2(tree);