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 c00c8771b86..07268ec93ed 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 @@ -1371,7 +1371,17 @@ 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_helper VPM_helper; + typedef typename VPM_helper::VPM SOTM_vpm2; + typedef typename VPM_helper::Tree_type Tree_type; + + 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)) { @@ -1382,16 +1392,20 @@ 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 (coplanar_patches_of_tm1.test(patch_id)) @@ -1401,12 +1415,13 @@ public: } 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); @@ -1414,9 +1429,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); @@ -1433,7 +1446,18 @@ 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_helper VPM_helper; + typedef typename VPM_helper::VPM SOTM_vpm1; + typedef typename VPM_helper::Tree_type Tree_type; + + 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)) { const std::size_t f_id = get(fids2, f); @@ -1443,14 +1467,17 @@ 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) @@ -1462,11 +1489,11 @@ public: } 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); @@ -1474,9 +1501,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 6c85bd8d8da..59c1e238972 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,207 @@ void copy_edge_mark(G&, No_mark&) {} // nothing to do + +// For exact side_of_triangle_mesh +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; +}; + +// For exact side_of_triangle_mesh +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_helper +{ + 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 + VPM get_vpm(const Node_id_map& node_ids, + const VertexPointMap& vpm, + const NodeVector& 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_helper +{ + typedef VertexPointMap VPM; + static + 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 // of input vertex point maps. template