diff --git a/BGL/include/CGAL/boost/graph/named_params_helper.h b/BGL/include/CGAL/boost/graph/named_params_helper.h index 4f589134439..8265365bb48 100644 --- a/BGL/include/CGAL/boost/graph/named_params_helper.h +++ b/BGL/include/CGAL/boost/graph/named_params_helper.h @@ -27,6 +27,20 @@ namespace CGAL { + namespace parameters + { + template + struct Is_default + { + typedef typename internal_np::Lookup_named_param_def < + Parameter, + NamedParameters, + internal_np::Param_not_found > ::type NP_type; + static const bool value = boost::is_same::value; + typedef CGAL::Boolean_tag type; + }; + } // end of parameters namespace + // forward declarations to avoid dependency to Solver_interface template class Default_diagonalize_traits; diff --git a/BGL/include/CGAL/boost/graph/parameters_interface.h b/BGL/include/CGAL/boost/graph/parameters_interface.h index 3e25ee8ff6c..14bdb1bdbca 100644 --- a/BGL/include/CGAL/boost/graph/parameters_interface.h +++ b/BGL/include/CGAL/boost/graph/parameters_interface.h @@ -102,6 +102,7 @@ CGAL_add_named_parameter(halfedges_keeper_t, halfedges_keeper, halfedges_keeper) CGAL_add_named_parameter(volume_threshold_t, volume_threshold, volume_threshold) CGAL_add_named_parameter(dry_run_t, dry_run, dry_run) CGAL_add_named_parameter(do_not_modify_t, do_not_modify, do_not_modify) +CGAL_add_named_parameter(non_manifold_feature_map_t, non_manifold_feature_map, non_manifold_feature_map) // List of named parameters that we use in the package 'Surface Mesh Simplification' CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost) diff --git a/BGL/test/BGL/test_cgal_bgl_named_params.cpp b/BGL/test/BGL/test_cgal_bgl_named_params.cpp index 1a2c89555b7..887a9098c13 100644 --- a/BGL/test/BGL/test_cgal_bgl_named_params.cpp +++ b/BGL/test/BGL/test_cgal_bgl_named_params.cpp @@ -99,6 +99,7 @@ void test(const NamedParameters& np) assert(get_parameter(np, CGAL::internal_np::do_simplify_border).v == 64); assert(get_parameter(np, CGAL::internal_np::do_not_modify).v == 65); assert(get_parameter(np, CGAL::internal_np::maximum_number_of_faces).v == 78910); + assert(get_parameter(np, CGAL::internal_np::non_manifold_feature_map).v == 60); // Named parameters that we use in the package 'Surface Mesh Simplification' assert(get_parameter(np, CGAL::internal_np::get_cost_policy).v == 34); @@ -207,6 +208,7 @@ void test(const NamedParameters& np) check_same_type<62>(get_parameter(np, CGAL::internal_np::halfedges_keeper)); check_same_type<64>(get_parameter(np, CGAL::internal_np::do_simplify_border)); check_same_type<78910>(get_parameter(np, CGAL::internal_np::maximum_number_of_faces)); + check_same_type<60>(get_parameter(np, CGAL::internal_np::non_manifold_feature_map)); // Named parameters that we use in the package 'Surface Mesh Simplification' check_same_type<34>(get_parameter(np, CGAL::internal_np::get_cost_policy)); @@ -325,6 +327,7 @@ int main() .throw_on_self_intersection(A<43>(43)) .clip_volume(A<44>(44)) .use_compact_clipper(A<45>(45)) + .non_manifold_feature_map(A<60>(60)) .apply_per_connected_component(A<46>(46)) .output_iterator(A<47>(47)) .erase_all_duplicates(A<48>(48)) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index ec8bb41756b..f3896b61eb1 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -84,6 +84,7 @@ create_single_source_cgal_program( "corefinement_mesh_union.cpp" ) create_single_source_cgal_program( "corefinement_mesh_union_and_intersection.cpp" ) create_single_source_cgal_program( "corefinement_mesh_union_with_attributes.cpp" ) create_single_source_cgal_program( "corefinement_polyhedron_union.cpp" ) +create_single_source_cgal_program( "corefine_non_manifold.cpp" ) create_single_source_cgal_program( "random_perturbation_SM_example.cpp" ) create_single_source_cgal_program( "corefinement_LCC.cpp") create_single_source_cgal_program( "detect_features_example.cpp" ) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefine_non_manifold.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefine_non_manifold.cpp new file mode 100644 index 00000000000..c3235d9f5f5 --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefine_non_manifold.cpp @@ -0,0 +1,232 @@ +#include +#include +#include +#include +#include + +#include + +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef K::Point_3 P; +typedef CGAL::Surface_mesh

Surface_mesh; + +namespace PMP = CGAL::Polygon_mesh_processing; + +int main() +{ +// polyline intersection with a non-manifold edge +{ + std::cout << "running polyline test\n"; + Surface_mesh tm1; + CGAL::make_quad(P(0,0,0), P(4,0,0), P(4,4,0), P(0,4,0), tm1); + CGAL::make_quad(P(0,-4,0), P(4,-4,0), P(4,0,0), P(0,0,0), tm1); + PMP::stitch_borders(tm1); + CGAL::make_quad(P(0,0,0), P(4,0,0), P(4,0,4), P(0,0,4), tm1); + CGAL::make_quad(P(0,0,0), P(4,0,0), P(4,0,-4), P(0,0,-4), tm1); + PMP::triangulate_faces(tm1); + + Surface_mesh tm2; + CGAL::make_quad(P(2,3,-2), P(2,-2,-2), P(2,-2,2), P(2,3,2), tm2); + PMP::triangulate_faces(tm2); + + PMP::Non_manifold_feature_map nm_map_1(tm1, get(boost::vertex_point, tm1)); + + std::vector< std::vector

> polylines; + PMP::surface_intersection(tm1, tm2, std::back_inserter(polylines), CGAL::parameters::non_manifold_feature_map(nm_map_1)); + + //dump polylines + std::ofstream output("intersection_polylines.cgal"); + for(const std::vector

& polyline : polylines) + { + output << polyline.size() << " "; + std::copy(polyline.begin(), polyline.end(),std::ostream_iterator

(output," ")); + output << "\n"; + } +} + +// simple case with only one non-manifold edge +{ + std::cout << "running corefinement test 1\n"; + Surface_mesh tm1; + CGAL::make_quad(P(0,0,0), P(4,0,0), P(4,4,0), P(0,4,0), tm1); + CGAL::make_quad(P(0,-4,0), P(4,-4,0), P(4,0,0), P(0,0,0), tm1); + PMP::stitch_borders(tm1); + CGAL::make_quad(P(0,0,0), P(4,0,0), P(4,0,4), P(0,0,4), tm1); + CGAL::make_quad(P(0,0,0), P(4,0,0), P(4,0,-4), P(0,0,-4), tm1); + PMP::triangulate_faces(tm1); + + Surface_mesh tm2; + CGAL::make_quad(P(2,3,-2), P(2,-2,-2), P(2,-2,2), P(2,3,2), tm2); + PMP::triangulate_faces(tm2); + + PMP::Non_manifold_feature_map nm_map_1(tm1, get(boost::vertex_point, tm1)); + + PMP::corefine(tm1, tm2, CGAL::parameters::non_manifold_feature_map(nm_map_1)); + + std::ofstream("t1_tm1_corefined.off") << tm1; + std::ofstream("t1_tm2_corefined.off") << tm2; +} + +// edge-edge intersection on a non-manifold edge +{ + std::cout << "running corefinement test 2\n"; + Surface_mesh tm1; + CGAL::make_quad(P(0,0,0), P(4,0,0), P(4,4,0), P(0,4,0), tm1); + CGAL::make_quad(P(0,-4,0), P(4,-4,0), P(4,0,0), P(0,0,0), tm1); + PMP::stitch_borders(tm1); + CGAL::make_quad(P(0,0,0), P(4,0,0), P(4,0,4), P(0,0,4), tm1); + CGAL::make_quad(P(0,0,0), P(4,0,0), P(4,0,-4), P(0,0,-4), tm1); + PMP::triangulate_faces(tm1); + + Surface_mesh tm2; + CGAL::make_quad(P(2,2,-2), P(2,-2,-2), P(2,-2,2), P(2,2,2), tm2); //TODO: test me + also test splitting the diagonal + PMP::triangulate_faces(tm2); + + PMP::Non_manifold_feature_map nm_map_1(tm1, get(boost::vertex_point, tm1)); + + PMP::corefine(tm1, tm2, CGAL::parameters::non_manifold_feature_map(nm_map_1)); + + std::ofstream("t2_tm1_corefined.off") << tm1; + std::ofstream("t2_tm2_corefined.off") << tm2; +} + +// coplanar edge and non-manifold edge +{ + std::cout << "running corefinement test 3\n"; + Surface_mesh tm1; + CGAL::make_quad(P(0,0,0), P(8,0,0), P(8,8,0), P(0,8,0), tm1); + CGAL::make_quad(P(8,0,0), P(12,0,0), P(12,8,0), P(8,8,0), tm1); + PMP::stitch_borders(tm1); + CGAL::make_quad(P(8,0,0), P(8,8,0), P(8,8,4), P(8,0,4), tm1); + CGAL::make_quad(P(8,0,0), P(8,8,0), P(8,8,-4), P(8,0,-4), tm1); + PMP::triangulate_faces(tm1); + + Surface_mesh tm2; + CGAL::make_quad(P(6,1,0), P(14,3,0), P(13,7,4), P(5,5,0), tm2); + CGAL::make_quad(P(6,1,0), P(5,5,0), P(5,5,2), P(6,1,3), tm2); + PMP::stitch_borders(tm2); + PMP::triangulate_faces(tm2); + + std::ofstream("t3_tm1.off") << tm1; + std::ofstream("t3_tm2.off") << tm2; + + PMP::Non_manifold_feature_map nm_map_1(tm1, get(boost::vertex_point, tm1)); + + PMP::corefine(tm1, tm2, CGAL::parameters::non_manifold_feature_map(nm_map_1)); + + std::ofstream("t3_tm1_corefined.off") << tm1; + std::ofstream("t3_tm2_corefined.off") << tm2; +} + +//TODO: add more tests nm-edge on vertex, nm-edge vs nm-edge, ... + + +// coplanar face and non-manifold edge +{ + std::cout << "running corefinement test 4\n"; + Surface_mesh tm1; + CGAL::make_quad(P(0,0,0), P(8,0,0), P(8,8,0), P(0,8,0), tm1); + CGAL::make_quad(P(8,0,0), P(12,0,0), P(12,8,0), P(8,8,0), tm1); + PMP::stitch_borders(tm1); + CGAL::make_quad(P(8,0,0), P(8,8,0), P(8,8,4), P(8,0,4), tm1); + CGAL::make_quad(P(8,0,0), P(8,8,0), P(8,8,-4), P(8,0,-4), tm1); + PMP::triangulate_faces(tm1); + + Surface_mesh tm2; + CGAL::make_quad(P(6,1,0), P(14,3,0), P(13,7,0), P(5,5,0), tm2); + CGAL::make_quad(P(6,1,0), P(5,5,0), P(5,5,2), P(6,1,3), tm2); + PMP::stitch_borders(tm2); + PMP::triangulate_faces(tm2); + + std::ofstream("t4_tm1.off") << tm1; + std::ofstream("t4_tm2.off") << tm2; + + PMP::Non_manifold_feature_map nm_map_1(tm1, get(boost::vertex_point, tm1)); + + PMP::corefine(tm1, tm2, CGAL::parameters::non_manifold_feature_map(nm_map_1)); + + std::ofstream("t4_tm1_corefined.off") << tm1; + std::ofstream("t4_tm2_corefined.off") << tm2; +} +// +// // coplanar face and non-manifold edge and regular intersection with incident face +{ + std::cout << "running corefinement test 5\n"; + Surface_mesh tm1; + CGAL::make_quad(P(0,0,0), P(4,0,0), P(4,8,0), P(0,8,0), tm1); + CGAL::make_quad(P(4,0,0), P(12,0,0), P(12,8,0), P(4,8,0), tm1); + PMP::stitch_borders(tm1); + CGAL::make_quad(P(4,0,0), P(4,8,0), P(4,8,4), P(4,0,4), tm1); + CGAL::make_quad(P(4,0,0), P(4,8,0), P(4,8,-4), P(4,0,-4), tm1); + PMP::triangulate_faces(tm1); + + Surface_mesh tm2; + CGAL::make_quad(P(2,1,0), P(14,3,0), P(13,7,0), P(5,5,0), tm2); + CGAL::make_quad(P(2,1,0), P(5,5,0), P(5,5,2), P(2,1,3), tm2); + PMP::stitch_borders(tm2); + PMP::triangulate_faces(tm2); + + std::ofstream("t5_tm1.off") << tm1; + std::ofstream("t5_tm2.off") << tm2; + + PMP::Non_manifold_feature_map nm_map_1(tm1, get(boost::vertex_point, tm1)); + + PMP::corefine(tm1, tm2, CGAL::parameters::non_manifold_feature_map(nm_map_1)); + + std::ofstream("t5_tm1_corefined.off") << tm1; + std::ofstream("t5_tm2_corefined.off") << tm2; +} + +// coplanar face and 2 non-manifold edges +{ + std::cout << "running corefinement test 6\n"; + Surface_mesh tm1; + CGAL::make_quad(P(0,0,0), P(8,0,0), P(8,8,0), P(0,8,0), tm1); + CGAL::make_quad(P(8,0,0), P(12,0,0), P(12,8,0), P(8,8,0), tm1); + PMP::stitch_borders(tm1); + CGAL::make_quad(P(8,0,0), P(8,8,0), P(8,8,4), P(8,0,4), tm1); + CGAL::make_quad(P(8,0,0), P(8,8,0), P(8,8,-4), P(8,0,-4), tm1); + PMP::triangulate_faces(tm1); + + Surface_mesh tm2; + CGAL::make_quad(P(6,1,0), P(14,3,0), P(13,7,0), P(5,5,0), tm2); + CGAL::make_quad(P(6,1,0), P(5,5,0), P(5,5,2), P(6,1,3), tm2); + PMP::stitch_borders(tm2); + CGAL::make_quad(P(6,1,0), P(5,5,0), P(5,5,-2), P(6,1,-3), tm2); + CGAL::make_quad(P(6,1,0), P(5,5,0), P(3,5,0), P(3,1,0), tm2); + PMP::triangulate_faces(tm2); + + std::ofstream("t6_tm1.off") << tm1; + std::ofstream("t6_tm2.off") << tm2; + + PMP::Non_manifold_feature_map nm_map_1(tm1, get(boost::vertex_point, tm1)); + PMP::Non_manifold_feature_map nm_map_2(tm2, get(boost::vertex_point, tm2)); + +#if 0 + std::vector< std::vector

> polylines; + PMP::surface_intersection(tm1, tm2, std::back_inserter(polylines), + CGAL::parameters::non_manifold_feature_map(nm_map_1), + CGAL::parameters::non_manifold_feature_map(nm_map_2)); + + //dump polylines + std::ofstream output("intersection_polylines.cgal"); + for(const std::vector

& polyline : polylines) + { + output << polyline.size() << " "; + std::copy(polyline.begin(), polyline.end(),std::ostream_iterator

(output," ")); + output << "\n"; + } +#else + PMP::corefine(tm1, tm2, CGAL::parameters::non_manifold_feature_map(nm_map_1), + CGAL::parameters::non_manifold_feature_map(nm_map_2)); + + std::ofstream("t6_tm1_corefined.off") << tm1; + std::ofstream("t6_tm2_corefined.off") << tm2; +#endif +} + + +} diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Non_manifold_feature_map.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Non_manifold_feature_map.h new file mode 100644 index 00000000000..419e330011e --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Non_manifold_feature_map.h @@ -0,0 +1,104 @@ +// Copyright (c) 2019 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Sebastien Loriot + + +#ifndef CGAL_POLYGON_MESH_PROCESSING_NON_MANIFOLD_FEATURE_MAP +#define CGAL_POLYGON_MESH_PROCESSING_NON_MANIFOLD_FEATURE_MAP + +namespace CGAL { +namespace Polygon_mesh_processing { + +//TODO: right now the base name parameter mechanism will make a deep copy, we probably want to avoid that +template +struct Non_manifold_feature_map +{ + typedef boost::graph_traits Graph_traits; + typedef typename Graph_traits::vertex_descriptor vertex_descriptor; + typedef typename Graph_traits::edge_descriptor edge_descriptor; + typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor; + typedef dynamic_edge_property_t Edge_to_id_tag; + typedef dynamic_vertex_property_t Vertex_to_id_tag; + typedef typename boost::property_map::type Edge_to_nm_id; + typedef typename boost::property_map::type Vertex_to_nm_id; + Edge_to_nm_id e_nm_id; + Vertex_to_nm_id v_nm_id; + std::vector< std::vector > non_manifold_edges; + std::vector< std::vector > non_manifold_vertices; + + Non_manifold_feature_map() + {} + + template + Non_manifold_feature_map(PolygonMesh& pm, Vpm vpm) + : e_nm_id(get(Edge_to_id_tag(), pm)) + , v_nm_id(get(Vertex_to_id_tag(), pm)) + { + typedef typename boost::property_traits::value_type Point_3; + + // detect non-manifold vertices + std::map > vertex_map; + for(vertex_descriptor vd : vertices(pm)) + { + put(v_nm_id, vd, std::size_t(-1)); // init map + vertex_map[get(vpm, vd)].push_back(vd); + } + + for(std::pair< const Point_3, + std::vector >& p : vertex_map) + { + if (p.second.size()!=1) + { + for (vertex_descriptor vd : p.second) + put(v_nm_id, vd, non_manifold_vertices.size()); + non_manifold_vertices.resize(non_manifold_vertices.size()+1); + // we steal the vertor from the map + p.second.swap(non_manifold_vertices.back()); + } + } + + // detect non-manifold edges + std::map< std::pair, std::vector > edge_map; + for(edge_descriptor ed : edges(pm)) + { + put(e_nm_id, ed, std::size_t(-1)); // init map + + halfedge_descriptor hd = halfedge(ed, pm); + + // an edge can be non-manifold only if both its vertices are non-manifold + if ( get(v_nm_id, source(hd, pm))==std::size_t(-1) || + get(v_nm_id, target(hd, pm))==std::size_t(-1) ) continue; + + const Point_3& src = get(vpm, source(ed, pm)); + const Point_3& tgt = get(vpm, target(ed, pm)); + // TODO: what to do with null edges? + if (src > tgt) + hd = opposite(hd, pm); + edge_map[ make_sorted_pair(src, tgt) ].push_back(edge(hd,pm)); + } + + for(std::pair< const std::pair, + std::vector >& p : edge_map) + { + if (p.second.size()!=1) + { + for (edge_descriptor ed : p.second) + put(e_nm_id, ed, non_manifold_edges.size()); + non_manifold_edges.resize(non_manifold_edges.size()+1); + // we steal the vertor from the map + p.second.swap(non_manifold_edges.back()); + } + } + } +}; + +} } // end of CGAL::Polygon_mesh_processing + +#endif diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/corefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/corefinement.h index e34c951e4ea..14a84753f10 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/corefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/corefinement.h @@ -779,15 +779,24 @@ corefine( TriangleMesh& tm1, > ::type User_visitor; User_visitor uv(choose_parameter(get_parameter(np1, internal_np::visitor))); + static const bool handle_non_manifold_features = + !parameters::Is_default::value || + !parameters::Is_default::value; + // surface intersection algorithm call typedef Corefinement::No_extra_output_from_corefinement Ob; typedef Corefinement::Surface_intersection_visitor_for_corefinement< - TriangleMesh, VPM1, VPM2, Ob, Ecm, User_visitor> Algo_visitor; + TriangleMesh, VPM1, VPM2, Ob, Ecm, User_visitor, false, handle_non_manifold_features> Algo_visitor; Ob ob; Ecm ecm(tm1,tm2,ecm1,ecm2); Corefinement::Intersection_of_triangle_meshes functor(tm1, tm2, vpm1, vpm2, Algo_visitor(uv,ob,ecm,const_mesh_ptr)); + + // Fill non-manifold feature maps if provided + functor.set_non_manifold_feature_map_1(parameters::get_parameter(np1, internal_np::non_manifold_feature_map)); + functor.set_non_manifold_feature_map_2(parameters::get_parameter(np2, internal_np::non_manifold_feature_map)); + functor(CGAL::Emptyset_iterator(), throw_on_self_intersection, true); } 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..01cdaa446b3 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 @@ -873,7 +873,7 @@ public: vpm1, vpm2, nodes) ) //p1==q2 { - CGAL_assertion( index_p1!=index_p2 || index_p1==Node_id((std::numeric_limits::max)()) ); + CGAL_assertion( index_p1!=index_p2 || index_p1==NID ); coplanar_patches_of_tm1.set(patch_id_p1); coplanar_patches_of_tm2.set(patch_id_q2); bool q1_is_between_p1p2 = sorted_around_edge( @@ -943,17 +943,17 @@ public: #endif //CGAL_COREFINEMENT_POLYHEDRA_DEBUG CGAL_assertion( - ( index_p1 == Node_id((std::numeric_limits::max)()) ? nodes.to_exact(get(vpm1,p1)): nodes.exact_node(index_p1) ) != - ( index_q1 == Node_id((std::numeric_limits::max)()) ? nodes.to_exact(get(vpm2,q1)): nodes.exact_node(index_q1) ) + ( index_p1 == NID ? nodes.to_exact(get(vpm1,p1)): nodes.exact_node(index_p1) ) != + ( index_q1 == NID ? nodes.to_exact(get(vpm2,q1)): nodes.exact_node(index_q1) ) && - ( index_p2 == Node_id((std::numeric_limits::max)()) ? nodes.to_exact(get(vpm1,p2)): nodes.exact_node(index_p2) ) != - ( index_q1 == Node_id((std::numeric_limits::max)()) ? nodes.to_exact(get(vpm2,q1)): nodes.exact_node(index_q1) ) + ( index_p2 == NID ? nodes.to_exact(get(vpm1,p2)): nodes.exact_node(index_p2) ) != + ( index_q1 == NID ? nodes.to_exact(get(vpm2,q1)): nodes.exact_node(index_q1) ) && - ( index_p1 == Node_id((std::numeric_limits::max)()) ? nodes.to_exact(get(vpm1,p1)): nodes.exact_node(index_p1) ) != - ( index_q2 == Node_id((std::numeric_limits::max)()) ? nodes.to_exact(get(vpm2,q2)): nodes.exact_node(index_q2) ) + ( index_p1 == NID ? nodes.to_exact(get(vpm1,p1)): nodes.exact_node(index_p1) ) != + ( index_q2 == NID ? nodes.to_exact(get(vpm2,q2)): nodes.exact_node(index_q2) ) && - ( index_p2 == Node_id((std::numeric_limits::max)()) ? nodes.to_exact(get(vpm1,p2)): nodes.exact_node(index_p2) ) != - ( index_q2 == Node_id((std::numeric_limits::max)()) ? nodes.to_exact(get(vpm2,q2)): nodes.exact_node(index_q2) ) + ( index_p2 == NID ? nodes.to_exact(get(vpm1,p2)): nodes.exact_node(index_p2) ) != + ( index_q2 == NID ? nodes.to_exact(get(vpm2,q2)): nodes.exact_node(index_q2) ) ); bool q1_is_between_p1p2 = sorted_around_edge( diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h index ca0d31915bb..a367f1eee99 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h @@ -107,7 +107,8 @@ template< class TriangleMesh, class OutputBuilder_ = Default, class EdgeMarkMapBind_ = Default, class UserVisitor_ = Default, - bool doing_autorefinement = false > + bool doing_autorefinement = false, + bool handle_non_manifold_features = false > class Surface_intersection_visitor_for_corefinement{ //default template parameters typedef typename Default::Get On_edge_map; //to keep the correspondance between node_id and vertex_handle in each mesh typedef std::vector Node_id_to_vertex; - typedef std::map Mesh_to_map_node; + typedef std::map Mesh_to_map_node; //to handle coplanar halfedge of polyhedra that are full in the intersection typedef std::multimap Node_to_target_of_hedge_map; typedef std::map Mesh_to_vertices_on_intersection_map; typedef boost::unordered_map Vertex_to_node_id; typedef std::map Mesh_to_vertex_to_node_id; + typedef Non_manifold_feature_map NM_features_map; // typedef for the CDT typedef Intersection_nodes INodes; @@ -153,7 +155,7 @@ private: typedef typename CDT::Vertex_handle CDT_Vertex_handle; // data members private: - // boost::dynamic_bitset<> non_manifold_nodes; + boost::dynamic_bitset<> is_node_on_boundary; // indicate if a vertex is a border vertex in tm1 or tm2 std::vector< std::vector > graph_of_constraints; boost::dynamic_bitset<> is_node_of_degree_one; //nb of intersection points between coplanar faces, see fixes XSL_TAG_CPL_VERT @@ -167,6 +169,9 @@ private: std::map< Node_id,std::set > coplanar_constraints; +// optional data members to handle non-manifold issues + std::map non_manifold_feature_maps; + //data members that require initialization in the constructor UserVisitor& user_visitor; OutputBuilder& output_builder; @@ -212,6 +217,70 @@ public: , const_mesh_ptr(const_mesh_ptr) {} + void + set_non_manifold_feature_map( + const TriangleMesh& tm, + const NM_features_map& nm) + { + non_manifold_feature_maps[&tm] = &nm; + } + + void copy_nodes_ids_for_non_manifold_features() + { + static const constexpr std::size_t NM_NID((std::numeric_limits::max)()); + + for(const std::pair& tm_and_nm : + non_manifold_feature_maps) + { + TriangleMesh* tm_ptr = const_cast(tm_and_nm.first); + // update nodes on edges + On_edge_map& on_edge_map = on_edge[tm_ptr]; + std::vector< std::pair > edges_to_copy; + for (const std::pair& ed_and_ids : on_edge_map) + { + std::size_t eid = get(tm_and_nm.second->e_nm_id, ed_and_ids.first); + if (eid!=NM_NID) + edges_to_copy.push_back(std::make_pair(eid,&(ed_and_ids.second))); + } + for(const std::pair& id_and_nodes : edges_to_copy) + { + const std::vector& nm_edges = + tm_and_nm.second->non_manifold_edges[id_and_nodes.first]; + CGAL_assertion( on_edge_map.count(nm_edges.front())==1 ); + + for (std::size_t i=1; i node_id + Vertex_to_node_id& vertex_to_node_id = mesh_to_vertex_to_node_id[tm_ptr]; + Node_to_target_of_hedge_map& vertices_on_inter = mesh_to_vertices_on_inter[tm_ptr]; + + std::vector< std::pair > vertices_to_add; + for (const typename std::pair& vd_and_id + : vertex_to_node_id) + { + std::size_t vid = get(tm_and_nm.second->v_nm_id, vd_and_id.first); + if (vid!=NM_NID) + vertices_to_add.push_back(std::make_pair(vd_and_id.first,vd_and_id.second)); + } + + for(const std::pair& vd_and_nid : vertices_to_add) + { + std::size_t vid = get(tm_and_nm.second->v_nm_id, vd_and_nid.first); + for(vertex_descriptor vd : tm_and_nm.second->non_manifold_vertices[vid]) + { + if (vd != vd_and_nid.first) + { + vertex_to_node_id.insert(std::make_pair(vd,vd_and_nid.second)); + output_builder.set_vertex_id(vd, vd_and_nid.second, *tm_ptr); + vertices_on_inter.insert(std::make_pair(vd_and_nid.second,halfedge(vd,*tm_ptr))); + } + } + } + } + } + template void annotate_graph(std::vector& graph) { @@ -220,14 +289,32 @@ public: is_node_of_degree_one.resize(nb_nodes); for(std::size_t node_id=0;node_id(&tm1); TriangleMesh* tm2_ptr = const_cast(&tm2); @@ -346,7 +426,7 @@ public: case ON_EDGE: //Edge intersected by an edge { on_edge[tm2_ptr][edge(h_2,tm2)].push_back(node_id); - // check_node_on_non_manifold_edge(node_id,h_2,tm2); + check_node_on_boundary_edge_case(node_id,h_2,tm2); } break; case ON_VERTEX: @@ -358,7 +438,7 @@ public: node_id_to_vertex.resize(node_id+1,Graph_traits::null_vertex()); node_id_to_vertex[node_id]=target(h_2,tm2); all_incident_faces_got_a_node_as_vertex(h_2,node_id,*tm2_ptr); - // check_node_on_non_manifold_vertex(node_id,h_2,tm2); + check_node_on_boundary_vertex_case(node_id,h_2,tm2); output_builder.set_vertex_id(target(h_2, tm2), node_id, tm2); } break; @@ -383,7 +463,7 @@ public: all_incident_faces_got_a_node_as_vertex(h_1,node_id, *tm1_ptr); // register the vertex in the output builder output_builder.set_vertex_id(target(h_1, tm1), node_id, tm1); - // check_node_on_non_manifold_vertex(node_id,h_1,tm1); + check_node_on_boundary_vertex_case(node_id,h_1,tm1); } else{ if ( is_source_coplanar ){ @@ -397,12 +477,34 @@ public: all_incident_faces_got_a_node_as_vertex(h_1_opp,node_id, *tm1_ptr); // register the vertex in the output builder output_builder.set_vertex_id(source(h_1, tm1), node_id, tm1); - // check_node_on_non_manifold_vertex(node_id,h_1_opp,tm1); + check_node_on_boundary_vertex_case(node_id,h_1_opp,tm1); } else{ //handle intersection on principal edge + typename std::map::iterator it_find = + non_manifold_feature_maps.find(&tm1); + if ( it_find != non_manifold_feature_maps.end() ) + { + // update h_1 if it is not the canonical non-manifold edge + // This is important to make sure intersection points on non-manifold + // edges are all connected for the same edge so that the redistribution + // on other edges does not overwrite some nodes. + // This update might be required in case of EDGE-EDGE intersection or + // COPLANAR intersection. + const NM_features_map& nm_features_map_1 = *it_find->second; + std::size_t eid1 = nm_features_map_1.non_manifold_edges.empty() + ? std::size_t(-1) + : get(nm_features_map_1.e_nm_id, edge(h_1, tm1)); + + if (eid1 != std::size_t(-1)) + { + if ( edge(h_1, tm1) != nm_features_map_1.non_manifold_edges[eid1].front() ) + h_1 = halfedge(nm_features_map_1.non_manifold_edges[eid1].front(), tm1); + } + } + on_edge[tm1_ptr][edge(h_1,tm1)].push_back(node_id); - // check_node_on_non_manifold_edge(node_id,h_1,tm1); + check_node_on_boundary_edge_case(node_id,h_1,tm1); } } } @@ -471,8 +573,9 @@ public: std::copy(begin,end,std::back_inserter(node_ids_array[it_id->second])); } - // Used by the autorefinement to re-set the id of nodes on the boundary of a - // face since another vertex (inside a face or on another edge) might have + // Used by the autorefinement and non-manifold edge handling to re-set + // the id of nodes on the boundary of a face since another vertex + // (inside a face or on another edge) might have // overwritten the vertex in node_id_to_vertex template void update_node_id_to_vertex_map(Node_id_to_vertex& node_id_to_vertex, @@ -564,7 +667,7 @@ public: // this condition ensures to consider only graph edges that are in // the same triangle if ( !points_on_triangle || it_vh!=id_to_CDT_vh.end() ){ - CGAL_assertion(doing_autorefinement || it_vh!=id_to_CDT_vh.end()); + CGAL_assertion(doing_autorefinement || handle_non_manifold_features || it_vh!=id_to_CDT_vh.end()); if (it_vh==id_to_CDT_vh.end()) continue; // needed for autorefinement (interior nodes) cdt.insert_constraint(vh,it_vh->second); constrained_edges.push_back(std::make_pair(id,id_n)); @@ -731,7 +834,7 @@ public: f_vertices[1]=it_fb->second.vertices[1]; f_vertices[2]=it_fb->second.vertices[2]; update_face_indices(f_vertices,f_indices,vertex_to_node_id); - if (doing_autorefinement) + if (doing_autorefinement || handle_non_manifold_features) it_fb->second.update_node_id_to_vertex_map(node_id_to_vertex, tm); } else{ @@ -777,7 +880,7 @@ public: { id_to_CDT_vh.insert( std::make_pair(f_indices[ik],triangle_vertices[ik])); - if (doing_autorefinement) + if (doing_autorefinement || handle_non_manifold_features) // update the current vertex in node_id_to_vertex // to match the one of the face node_id_to_vertex[f_indices[ik]]=f_vertices[ik]; @@ -952,6 +1055,8 @@ public: const VertexPointMap1& vpm1, const VertexPointMap2& vpm2) { + copy_nodes_ids_for_non_manifold_features(); + nodes.all_nodes_created(); TriangleMesh* tm1_ptr = const_cast(&tm1); @@ -984,15 +1089,19 @@ public: // Face_boundaries& face_boundaries=mesh_to_face_boundaries[&tm]; Node_to_target_of_hedge_map& nodes_to_hedge=it->second; + + // iterate on the vertices that are on the intersection between the input meshes for(typename Node_to_target_of_hedge_map::iterator it_node_2_hedge=nodes_to_hedge.begin(); it_node_2_hedge!=nodes_to_hedge.end(); ++it_node_2_hedge) { Node_id node_id_of_first=it_node_2_hedge->first; + // look for neighbors of the current node in the intersection graph std::vector& neighbors=graph_of_constraints[node_id_of_first]; if ( !neighbors.empty() ) { + // for all neighbors look for input vertices that are also on the intersection for(Node_id node_id : neighbors) { //if already done for the opposite @@ -1010,7 +1119,7 @@ public: target(it_node_2_hedge_two->second,tm) ) { hedge=opposite(next(hedge,tm),tm); - if (tm1_ptr==tm2_ptr && hedge==start) + if ((doing_autorefinement || handle_non_manifold_features) && hedge==start) { ++it_node_2_hedge_two; // we are using a multimap and // the halfedge we are looking for @@ -1066,7 +1175,7 @@ public: triangulate_intersected_faces(it, vpm2, nodes, mesh_to_face_boundaries); } - nodes.finalize(); + nodes.finalize(mesh_to_node_id_to_vertex); // additional operations output_builder(nodes, 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..26881b21f61 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 @@ -521,28 +521,38 @@ struct Patch_container{ typedef typename GT::face_descriptor face_descriptor; Patch_description& patch=this->operator[](i); - out << "OFF\n" << patch.interior_vertices.size() + - patch.shared_edges.size(); - out << " " << patch.faces.size() << " 0\n"; + + std::stringstream ss; std::map vertexid; int id=0; for(vertex_descriptor vh : patch.interior_vertices) { vertexid[vh]=id++; - out << get(vertex_point, pm, vh) << "\n"; + ss << get(vertex_point, pm, vh) << "\n"; } for(halfedge_descriptor hh : patch.shared_edges) { - vertexid[target(hh, pm)]=id++; - out << get(vertex_point, pm, target(hh, pm)) << "\n"; + if( vertexid.insert( std::make_pair(target(hh,pm), id) ).second ) + { + ss << get(vertex_point, pm, target(hh, pm)) << "\n"; + ++id; + } + if( vertexid.insert( std::make_pair(source(hh,pm), id) ).second ) + { + ss << get(vertex_point, pm, source(hh, pm)) << "\n"; + ++id; + } } + out << "OFF\n" << id << " " << patch.faces.size() << " 0\n"; + out << ss.str(); + for(face_descriptor f : patch.faces) { - out << "3 " << vertexid[source(halfedge(f,pm),pm)] << - " " << vertexid[target(halfedge(f,pm),pm)] << - " " << vertexid[target(next(halfedge(f,pm),pm),pm)] << "\n"; + out << "3 " << vertexid.at(source(halfedge(f,pm),pm)) << + " " << vertexid.at(target(halfedge(f,pm),pm)) << + " " << vertexid.at(target(next(halfedge(f,pm),pm),pm)) << "\n"; } return out; @@ -559,6 +569,16 @@ struct Patch_container{ dump_patch(i, output); } } + void dump_patches(std::string prefix) + { + for (std::size_t i=0;i #include #include +#include #include #include @@ -98,9 +99,12 @@ struct Default_surface_intersection_visitor{ // If we implement a predicate only test, we can get rid of it. static const bool Predicates_on_constructions_needed = doing_autorefinement; static const bool do_need_vertex_graph = false; + void set_non_manifold_feature_map( + const TriangleMesh&, + const Non_manifold_feature_map&) + {} }; - struct Node_id_set { typedef std::size_t Node_id; @@ -187,6 +191,9 @@ class Intersection_of_triangle_meshes Node_visitor visitor; Faces_to_nodes_map f_to_node; //Associate a pair of triangles to their intersection points std::vector extra_terminal_nodes; //used only for autorefinement + Non_manifold_feature_map non_manifold_feature_map_1, + non_manifold_feature_map_2; + static const constexpr std::size_t NM_NID = (std::numeric_limits::max)(); CGAL_assertion_code(bool doing_autorefinement;) // member functions @@ -195,6 +202,7 @@ class Intersection_of_triangle_meshes const TriangleMesh& tm_e, const VPMF& vpm_f, const VPME& vpm_e, + const Non_manifold_feature_map& non_manifold_feature_map, bool throw_on_self_intersection) { std::vector face_boxes, edge_boxes; @@ -215,15 +223,38 @@ class Intersection_of_triangle_meshes edge_boxes.reserve(num_edges(tm_e)); edge_boxes_ptr.reserve(num_edges(tm_e)); - for(edge_descriptor ed : edges(tm_e)) - { - halfedge_descriptor h=halfedge(ed,tm_e); - edge_boxes.push_back( Box( - get(vpm_e,source(h,tm_e)).bbox() + - get(vpm_e,target(h,tm_e)).bbox(), - h ) ); - edge_boxes_ptr.push_back( &edge_boxes.back() ); - } + if (non_manifold_feature_map.non_manifold_edges.empty()) + // general manifold case + for(edge_descriptor ed : edges(tm_e)) + { + halfedge_descriptor h=halfedge(ed,tm_e); + edge_boxes.push_back( Box( + get(vpm_e,source(h,tm_e)).bbox() + + get(vpm_e,target(h,tm_e)).bbox(), + h ) ); + edge_boxes_ptr.push_back( &edge_boxes.back() ); + } + else + // non-manifold case + for(edge_descriptor ed : edges(tm_e)) + { + std::size_t eid=get(non_manifold_feature_map.e_nm_id, ed); + halfedge_descriptor h=halfedge(ed,tm_e); + // insert only one copy of a non-manifold edge + if (eid!=NM_NID) + { + if (non_manifold_feature_map.non_manifold_edges[eid].front()!=ed) + continue; + else + // make sure the halfedge used is consistant with stored one + h = halfedge(non_manifold_feature_map.non_manifold_edges[eid].front(), tm_e); + } + edge_boxes.push_back( Box( + get(vpm_e,source(h,tm_e)).bbox() + + get(vpm_e,target(h,tm_e)).bbox(), + h ) ); + edge_boxes_ptr.push_back( &edge_boxes.back() ); + } /// \todo experiments different cutoff values std::ptrdiff_t cutoff = 2 * std::ptrdiff_t( @@ -314,18 +345,34 @@ class Intersection_of_triangle_meshes get_or_create_node(const Cpl_inter_pt& ipt, Node_id& current_node, std::map& coplanar_node_map, - const TriangleMesh& tm1, - const TriangleMesh& tm2) + const Non_manifold_feature_map& nm_features_map_1, + const Non_manifold_feature_map& nm_features_map_2, + const TriangleMesh& tm1, + const TriangleMesh& tm2) { halfedge_descriptor h1=graph_traits::null_halfedge(),h2=h1; switch(ipt.type_1){ case ON_VERTEX: - h1=halfedge(target(ipt.info_1,tm1),tm1); + { + std::size_t vid1 = nm_features_map_1.non_manifold_vertices.empty() + ? NM_NID + : get(nm_features_map_1.v_nm_id, target(ipt.info_1,tm1)); + if (vid1==NM_NID) + h1=halfedge(target(ipt.info_1,tm1),tm1); + else + h1=halfedge(nm_features_map_1.non_manifold_vertices[vid1][0],tm1); + } break; case ON_EDGE : { - h1=opposite(ipt.info_1,tm1); - if (h1>ipt.info_1) h1=ipt.info_1; + std::size_t eid1 = nm_features_map_1.non_manifold_edges.empty() + ? NM_NID + : get(nm_features_map_1.e_nm_id, edge(ipt.info_1,tm1)); + if (eid1==NM_NID) + h1=ipt.info_1; + else + h1=halfedge(nm_features_map_1.non_manifold_edges[eid1][0],tm1); + h1=(std::max)(h1, opposite(h1, tm1)); } break; case ON_FACE : @@ -335,12 +382,26 @@ class Intersection_of_triangle_meshes } switch(ipt.type_2){ case ON_VERTEX: - h2=halfedge(target(ipt.info_2,tm2),tm2); + { + std::size_t vid2 = nm_features_map_2.non_manifold_vertices.empty() + ? NM_NID + : get(nm_features_map_2.v_nm_id, target(ipt.info_2,tm2)); + if (vid2==NM_NID) + h2=halfedge(target(ipt.info_2,tm2),tm2); + else + h2=halfedge(nm_features_map_2.non_manifold_vertices[vid2][0],tm2); + } break; case ON_EDGE : { - h2=opposite(ipt.info_2,tm2); - if (h2>ipt.info_2) h2=ipt.info_2; + std::size_t eid2 = nm_features_map_2.non_manifold_edges.empty() + ? NM_NID + : get(nm_features_map_2.e_nm_id, edge(ipt.info_2,tm2)); + if (eid2==NM_NID) + h2=ipt.info_2; + else + h2=halfedge(nm_features_map_2.non_manifold_edges[eid2][0],tm2); + h2=(std::max)(h2, opposite(h2, tm2)); } break; case ON_FACE : @@ -363,14 +424,14 @@ class Intersection_of_triangle_meshes } void add_intersection_point_to_face_and_all_edge_incident_faces(face_descriptor f_1, - halfedge_descriptor e_2, + halfedge_descriptor h_2, const TriangleMesh& tm1, const TriangleMesh& tm2, Node_id node_id) { - if (!is_border(e_2, tm2)) + if (!is_border(h_2, tm2)) { - face_descriptor f_2 = face(e_2, tm2); + face_descriptor f_2 = face(h_2, tm2); if(&tm1!=&tm2 || f_1!=f_2) { Face_pair face_pair = &tm1==&tm2 ? make_sorted_pair(f_1,f_2): @@ -381,10 +442,10 @@ class Intersection_of_triangle_meshes f_to_node[Face_pair_and_int(face_pair,0)].insert(node_id); } } - e_2 = opposite(e_2, tm2); - if (!is_border(e_2, tm2)) + h_2 = opposite(h_2, tm2); + if (!is_border(h_2, tm2)) { - face_descriptor f_2 = face(e_2, tm2); + face_descriptor f_2 = face(h_2, tm2); if(&tm1!=&tm2 || f_1!=f_2) { Face_pair face_pair = &tm1==&tm2 ? make_sorted_pair(f_1,f_2): @@ -465,6 +526,7 @@ class Intersection_of_triangle_meshes halfedge_descriptor f_2, const TriangleMesh& tm1, const TriangleMesh& tm2, + const Non_manifold_feature_map& nm_features_map_1, Node_id node_id, bool is_new_node) { @@ -475,42 +537,78 @@ class Intersection_of_triangle_meshes ? stm_edge_to_ltm_faces : ltm_edge_to_stm_faces; - for(halfedge_descriptor h_1 : - halfedges_around_target(v_1,tm1)) - { - add_intersection_point_to_face_and_all_edge_incident_faces(face(f_2,tm2),h_1,tm2,tm1,node_id); - typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(h_1,tm1)); - if (it_ets!=tm1_edge_to_tm2_faces.end()) it_ets->second.erase(face(f_2,tm2)); - } + std::vector tmp_vertices_1(1, target(v_1, tm1)); + + std::size_t vid1 = nm_features_map_1.non_manifold_vertices.empty() + ? NM_NID + : get(nm_features_map_1.v_nm_id, target(v_1, tm1)); + + const std::vector& vertices_1 = vid1==NM_NID + ? tmp_vertices_1 + : nm_features_map_1.non_manifold_vertices[vid1]; + for(vertex_descriptor v1 : vertices_1) + for(halfedge_descriptor h_1 : + halfedges_around_target(v1,tm1)) + { + add_intersection_point_to_face_and_all_edge_incident_faces(face(f_2,tm2),h_1,tm2,tm1,node_id); + typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(h_1,tm1)); + if (it_ets!=tm1_edge_to_tm2_faces.end()) it_ets->second.erase(face(f_2,tm2)); + } } void handle_coplanar_case_VERTEX_EDGE(halfedge_descriptor v_1, - halfedge_descriptor e_2, + halfedge_descriptor h_2, const TriangleMesh& tm1, const TriangleMesh& tm2, + const Non_manifold_feature_map& nm_features_map_1, + const Non_manifold_feature_map& nm_features_map_2, Node_id node_id, bool is_new_node) { if(is_new_node) - visitor.new_node_added(node_id,ON_VERTEX,e_2,v_1,tm2,tm1,false,false); + visitor.new_node_added(node_id,ON_VERTEX,h_2,v_1,tm2,tm1,false,false); Edge_to_faces& tm1_edge_to_tm2_faces = &tm1 <= &tm2 ? stm_edge_to_ltm_faces : ltm_edge_to_stm_faces; - for(halfedge_descriptor h_1 : - halfedges_around_target(v_1,tm1)) - { - typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(h_1,tm1)); - Face_set* fset = (it_ets!=tm1_edge_to_tm2_faces.end())?&(it_ets->second):nullptr; - cip_handle_case_edge(node_id,fset,h_1,e_2,tm1,tm2); - } + std::vector tmp_vertices_1(1, target(v_1, tm1)); + + std::size_t vid1 = nm_features_map_1.non_manifold_vertices.empty() + ? NM_NID + : get(nm_features_map_1.v_nm_id, target(v_1, tm1)); + + const std::vector& vertices_1 = vid1==NM_NID + ? tmp_vertices_1 + : nm_features_map_1.non_manifold_vertices[vid1]; + + std::vector tmp_edges_2(1, edge(h_2, tm2)); + + std::size_t eid2 = nm_features_map_2.non_manifold_edges.empty() + ? NM_NID + : get(nm_features_map_2.e_nm_id, edge(h_2, tm2)); + + const std::vector& edges_2 = eid2==NM_NID + ? tmp_edges_2 + : nm_features_map_2.non_manifold_edges[eid2]; + + for (vertex_descriptor v1 : vertices_1) + for(halfedge_descriptor h_1 : + halfedges_around_target(v1,tm1)) + { + typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(h_1,tm1)); + Face_set* fset = (it_ets!=tm1_edge_to_tm2_faces.end())?&(it_ets->second):nullptr; + for (edge_descriptor e2 : edges_2) + cip_handle_case_edge(node_id,fset,h_1,halfedge(e2, tm2),tm1,tm2); + } } void handle_coplanar_case_VERTEX_VERTEX(halfedge_descriptor v_1, halfedge_descriptor v_2, const TriangleMesh& tm1, const TriangleMesh& tm2, + const Non_manifold_feature_map& nm_features_map_1, + const Non_manifold_feature_map& nm_features_map_2, Node_id node_id, bool is_new_node) { @@ -521,13 +619,33 @@ class Intersection_of_triangle_meshes ? stm_edge_to_ltm_faces : ltm_edge_to_stm_faces; - for(halfedge_descriptor h_1 : - halfedges_around_target(v_1,tm1)) - { - typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(h_1,tm1)); - Face_set* fset = (it_ets!=tm1_edge_to_tm2_faces.end())?&(it_ets->second):nullptr; - cip_handle_case_vertex(node_id,fset,h_1,v_2,tm1,tm2); - } + std::vector tmp_vertices_1(1, target(v_1, tm1)), + tmp_vertices_2(1, target(v_2, tm2)); + + std::size_t vid1 = nm_features_map_1.non_manifold_vertices.empty() + ? NM_NID + : get(nm_features_map_1.v_nm_id, target(v_1, tm1)); + + std::size_t vid2 = nm_features_map_2.non_manifold_vertices.empty() + ? NM_NID + : get(nm_features_map_2.v_nm_id, target(v_2, tm2)); + + const std::vector& vertices_1 = vid1==NM_NID + ? tmp_vertices_1 + : nm_features_map_1.non_manifold_vertices[vid1]; + const std::vector& vertices_2 = vid2==NM_NID + ? tmp_vertices_2 + : nm_features_map_2.non_manifold_vertices[vid2]; + + for (vertex_descriptor v1 : vertices_1) + for(halfedge_descriptor h_1 : + halfedges_around_target(v1,tm1)) + { + typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(h_1,tm1)); + Face_set* fset = (it_ets!=tm1_edge_to_tm2_faces.end())?&(it_ets->second):nullptr; + for (vertex_descriptor v2 : vertices_2) + cip_handle_case_vertex(node_id,fset,h_1,halfedge(v2, tm2),tm1,tm2); + } } template @@ -535,7 +653,9 @@ class Intersection_of_triangle_meshes const TriangleMesh& tm1, const TriangleMesh& tm2, const VPM1& vpm1, - const VPM2& vpm2) + const VPM2& vpm2, + const Non_manifold_feature_map& nm_features_map_1, + const Non_manifold_feature_map& nm_features_map_2) { CGAL_assertion( &tm1 < &tm2 || &tm1==&tm2 ); @@ -569,7 +689,7 @@ class Intersection_of_triangle_meshes Node_id node_id; bool is_new_node; std::tie(node_id, is_new_node) = - get_or_create_node(ipt,current_node,coplanar_node_map,tm1,tm2); + get_or_create_node(ipt,current_node,coplanar_node_map,nm_features_map_1,nm_features_map_2,tm1,tm2); cpln_nodes.push_back(node_id); switch(ipt.type_1){ @@ -577,13 +697,13 @@ class Intersection_of_triangle_meshes { switch(ipt.type_2){ case ON_VERTEX: - handle_coplanar_case_VERTEX_VERTEX(ipt.info_1,ipt.info_2,tm1,tm2,node_id,is_new_node); + handle_coplanar_case_VERTEX_VERTEX(ipt.info_1,ipt.info_2,tm1,tm2,nm_features_map_1,nm_features_map_2,node_id,is_new_node); break; case ON_EDGE: - handle_coplanar_case_VERTEX_EDGE(ipt.info_1,ipt.info_2,tm1,tm2,node_id,is_new_node); + handle_coplanar_case_VERTEX_EDGE(ipt.info_1,ipt.info_2,tm1,tm2,nm_features_map_1,nm_features_map_2,node_id,is_new_node); break; case ON_FACE: - handle_coplanar_case_VERTEX_FACE(ipt.info_1,ipt.info_2,tm1,tm2,node_id,is_new_node); + handle_coplanar_case_VERTEX_FACE(ipt.info_1,ipt.info_2,tm1,tm2,nm_features_map_1,node_id,is_new_node); break; default: CGAL_error_msg("Should not get there!"); } @@ -593,15 +713,33 @@ class Intersection_of_triangle_meshes { switch(ipt.type_2){ case ON_VERTEX: - handle_coplanar_case_VERTEX_EDGE(ipt.info_2,ipt.info_1,tm2,tm1,node_id,is_new_node); + handle_coplanar_case_VERTEX_EDGE(ipt.info_2,ipt.info_1,tm2,tm1,nm_features_map_2,nm_features_map_1,node_id,is_new_node); break; case ON_EDGE: { + std::vector tmp_edges_1(1, edge(ipt.info_1,tm1)), + tmp_edges_2(1, edge(ipt.info_2,tm2)); + std::size_t eid1 = nm_features_map_1.non_manifold_edges.empty() + ? NM_NID + : get(nm_features_map_1.e_nm_id, edge(ipt.info_1, tm1)); + std::size_t eid2 = nm_features_map_2.non_manifold_edges.empty() + ? NM_NID + : get(nm_features_map_2.e_nm_id, edge(ipt.info_2, tm2)); + const std::vector& edges_1 = eid1==NM_NID + ? tmp_edges_1 + : nm_features_map_1.non_manifold_edges[eid1]; + const std::vector& edges_2 = eid2==NM_NID + ? tmp_edges_2 + : nm_features_map_2.non_manifold_edges[eid2]; if(is_new_node) - visitor.new_node_added(node_id,ON_EDGE,ipt.info_1,ipt.info_2,tm1,tm2,false,false); - typename Edge_to_faces::iterator it_ets=stm_edge_to_ltm_faces.find(edge(ipt.info_1,tm1)); - Face_set* fset = (it_ets!=stm_edge_to_ltm_faces.end())?&(it_ets->second):nullptr; - cip_handle_case_edge(node_id,fset,ipt.info_1,ipt.info_2,tm1,tm2); + visitor.new_node_added(node_id,ON_EDGE,halfedge(edges_1.front(), tm1),halfedge(edges_2.front(),tm2),tm1,tm2,false,false); + for(edge_descriptor e1 : edges_1) + for(edge_descriptor e2 : edges_2) + { + typename Edge_to_faces::iterator it_ets=stm_edge_to_ltm_faces.find(e1); + Face_set* fset = (it_ets!=stm_edge_to_ltm_faces.end())?&(it_ets->second):nullptr; + cip_handle_case_edge(node_id,fset,halfedge(e1,tm1),halfedge(e2,tm2),tm1,tm2); + } } break; default: CGAL_error_msg("Should not get there!"); @@ -611,7 +749,7 @@ class Intersection_of_triangle_meshes case ON_FACE: { CGAL_assertion(ipt.type_2==ON_VERTEX); - handle_coplanar_case_VERTEX_FACE(ipt.info_2,ipt.info_1,tm2,tm1,node_id,is_new_node); + handle_coplanar_case_VERTEX_FACE(ipt.info_2,ipt.info_1,tm2,tm1,nm_features_map_2,node_id,is_new_node); } break; default: CGAL_error_msg("Should not get there!"); @@ -667,6 +805,8 @@ class Intersection_of_triangle_meshes const TriangleMesh& tm2, const VPM1& vpm1, const VPM2& vpm2, + const Non_manifold_feature_map& nm_features_map_1, + const Non_manifold_feature_map& nm_features_map_2, Node_id& current_node) { typedef std::tuple Inter_type; @@ -676,6 +816,7 @@ class Intersection_of_triangle_meshes it!=tm1_edge_to_tm2_faces.end();++it) { edge_descriptor e_1=it->first; + halfedge_descriptor h_1=halfedge(e_1,tm1); Face_set& fset=it->second; while (!fset.empty()){ @@ -687,18 +828,93 @@ class Intersection_of_triangle_meshes //handle degenerate case: one extremity of edge belong to f_2 std::vector all_edges; if ( std::get<3>(res) ) // is edge target in triangle plane - std::copy(halfedges_around_target(h_1,tm1).first, - halfedges_around_target(h_1,tm1).second, - std::back_inserter(all_edges)); + { + if (!nm_features_map_1.non_manifold_edges.empty()) + { + std::size_t vid1 = get(nm_features_map_1.v_nm_id, target(h_1, tm1)); + if (vid1 != NM_NID) + { + for (vertex_descriptor vd : nm_features_map_1.non_manifold_vertices[vid1]) + { + std::copy(halfedges_around_target(vd,tm1).first, + halfedges_around_target(vd,tm1).second, + std::back_inserter(all_edges)); + } + if (all_edges.front()!=h_1) + { + // restore expected property + typename std::vector::iterator pos = + std::find(all_edges.begin(), all_edges.end(), h_1); + CGAL_assertion(pos!=all_edges.end()); + std::swap(*pos, all_edges.front()); + } + } + else + std::copy(halfedges_around_target(h_1,tm1).first, + halfedges_around_target(h_1,tm1).second, + std::back_inserter(all_edges)); + } + else + std::copy(halfedges_around_target(h_1,tm1).first, + halfedges_around_target(h_1,tm1).second, + std::back_inserter(all_edges)); + } else{ if ( std::get<2>(res) ) // is edge source in triangle plane + { + if (!nm_features_map_1.non_manifold_edges.empty()) + { + std::size_t vid1 = get(nm_features_map_1.v_nm_id, source(h_1, tm1)); + if (vid1 != NM_NID) + { + for (vertex_descriptor vd : nm_features_map_1.non_manifold_vertices[vid1]) + { + std::copy(halfedges_around_source(vd,tm1).first, + halfedges_around_source(vd,tm1).second, + std::back_inserter(all_edges)); + } + if (all_edges.front()!=h_1) + { + // restore expected property + typename std::vector::iterator pos = + std::find(all_edges.begin(), all_edges.end(), h_1); + CGAL_assertion(pos!=all_edges.end()); + std::swap(*pos, all_edges.front()); + } + } + else + std::copy(halfedges_around_source(h_1,tm1).first, + halfedges_around_source(h_1,tm1).second, + std::back_inserter(all_edges)); + } + else std::copy(halfedges_around_source(h_1,tm1).first, halfedges_around_source(h_1,tm1).second, std::back_inserter(all_edges)); + } else + { all_edges.push_back(h_1); + edge_descriptor e_1 = edge(h_1, tm1); + if (!nm_features_map_1.non_manifold_edges.empty()) + { + std::size_t eid1 = get(nm_features_map_1.e_nm_id, e_1); + if (eid1 != NM_NID) + { + CGAL_assertion( nm_features_map_1.non_manifold_edges[eid1][0]==e_1 ); + for (std::size_t k=1; + k(res); + + std::size_t eid2 = nm_features_map_2.non_manifold_edges.empty() + ? NM_NID + : get(nm_features_map_2.e_nm_id, edge(h_2, tm2)); + + if (eid2!=NM_NID) + h_2 = halfedge(nm_features_map_2.non_manifold_edges[eid2].front(), tm2); + visitor.new_node_added(node_id,ON_EDGE,h_1,h_2,tm1,tm2,std::get<3>(res),std::get<2>(res)); for (;it_edge!=all_edges.end();++it_edge){ if ( it_edge!=all_edges.begin() ){ typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(*it_edge,tm1)); Face_set* fset_bis = (it_ets!=tm1_edge_to_tm2_faces.end())?&(it_ets->second):nullptr; - cip_handle_case_edge(node_id,fset_bis,*it_edge,h_2,tm1,tm2); + if( eid2 == NM_NID ) + cip_handle_case_edge(node_id,fset_bis,*it_edge,h_2,tm1,tm2); + else + { + for (edge_descriptor e2 : nm_features_map_2.non_manifold_edges[eid2]) + cip_handle_case_edge(node_id,fset_bis,*it_edge,halfedge(e2, tm2),tm1,tm2); + } } else - cip_handle_case_edge(node_id,&fset,*it_edge,h_2,tm1,tm2); + { + if( eid2 == NM_NID ) + cip_handle_case_edge(node_id,&fset,*it_edge,h_2,tm1,tm2); + else + for (edge_descriptor e2 : nm_features_map_2.non_manifold_edges[eid2]) + cip_handle_case_edge(node_id,&fset,*it_edge,halfedge(e2, tm2),tm1,tm2); + } } } // end case ON_EDGE break; @@ -767,14 +1005,27 @@ class Intersection_of_triangle_meshes nodes.add_new_node(get(vpm2, target(h_2,tm2))); //we use the original vertex to create the node //before it was ON_FACE but do not remember why, probably a bug... visitor.new_node_added(node_id,ON_VERTEX,h_1,h_2,tm1,tm2,std::get<3>(res),std::get<2>(res)); + + std::size_t vid2 = nm_features_map_2.non_manifold_vertices.empty() + ? NM_NID + : get(nm_features_map_2.v_nm_id, target(h_2, tm2)); + for (;it_edge!=all_edges.end();++it_edge){ if ( it_edge!=all_edges.begin() ){ typename Edge_to_faces::iterator it_ets=tm1_edge_to_tm2_faces.find(edge(*it_edge,tm1)); Face_set* fset_bis = (it_ets!=tm1_edge_to_tm2_faces.end())?&(it_ets->second):nullptr; - cip_handle_case_vertex(node_id,fset_bis,*it_edge,h_2,tm1,tm2); + if( vid2 == NM_NID ) + cip_handle_case_vertex(node_id,fset_bis,*it_edge,h_2,tm1,tm2); + else + for (vertex_descriptor vd2 : nm_features_map_2.non_manifold_vertices[vid2]) + cip_handle_case_vertex(node_id,fset_bis,*it_edge,halfedge(vd2, tm2),tm1,tm2); } else - cip_handle_case_vertex(node_id,&fset,*it_edge,h_2,tm1,tm2); + if( vid2 == NM_NID ) + cip_handle_case_vertex(node_id,&fset,*it_edge,h_2,tm1,tm2); + else + for (vertex_descriptor vd2 : nm_features_map_2.non_manifold_vertices[vid2]) + cip_handle_case_vertex(node_id,&fset,*it_edge,halfedge(vd2, tm2),tm1,tm2); } } // end case ON_VERTEX break; @@ -1104,7 +1355,7 @@ class Intersection_of_triangle_meshes } else{ CGAL_assertion(segment.size()==1); - isolated_point_seen=true; + isolated_point_seen=true; // NOT TRUE CAN BE END POINT OF POLYLINE FALLING ONTO AN INPUT EDGE } } @@ -1290,6 +1541,20 @@ public: CGAL_assertion_code( doing_autorefinement=true; ) } +// setting maps of non manifold features + void set_non_manifold_feature_map_1(internal_np::Param_not_found){} + void set_non_manifold_feature_map_2(internal_np::Param_not_found){} + void set_non_manifold_feature_map_1(const Non_manifold_feature_map& m) + { + non_manifold_feature_map_1=m; + visitor.set_non_manifold_feature_map(nodes.tm1, non_manifold_feature_map_1); + } + void set_non_manifold_feature_map_2(const Non_manifold_feature_map& m) + { + non_manifold_feature_map_2=m; + visitor.set_non_manifold_feature_map(nodes.tm2, non_manifold_feature_map_2); + } + template OutputIterator operator()(OutputIterator output, bool throw_on_self_intersection, @@ -1302,18 +1567,18 @@ public: const VertexPointMap1& vpm1=nodes.vpm1; const VertexPointMap2& vpm2=nodes.vpm2; - filter_intersections(tm1, tm2, vpm1, vpm2, throw_on_self_intersection); - filter_intersections(tm2, tm1, vpm2, vpm1, throw_on_self_intersection); + filter_intersections(tm1, tm2, vpm1, vpm2, non_manifold_feature_map_2, throw_on_self_intersection); + filter_intersections(tm2, tm1, vpm2, vpm1, non_manifold_feature_map_1, throw_on_self_intersection); Node_id current_node((std::numeric_limits::max)()); CGAL_assertion(current_node+1==0); - +// TODO: handle non-manifold edges in coplanar #ifndef DO_NOT_HANDLE_COPLANAR_FACES //first handle coplanar triangles if (&tm1<&tm2) - compute_intersection_of_coplanar_faces(current_node, tm1, tm2, vpm1, vpm2); + compute_intersection_of_coplanar_faces(current_node, tm1, tm2, vpm1, vpm2, non_manifold_feature_map_1, non_manifold_feature_map_2); else - compute_intersection_of_coplanar_faces(current_node, tm2, tm1, vpm2, vpm1); + compute_intersection_of_coplanar_faces(current_node, tm2, tm1, vpm2, vpm1, non_manifold_feature_map_2, non_manifold_feature_map_1); visitor.set_number_of_intersection_points_from_coplanar_faces(current_node+1); if (!coplanar_faces.empty()) @@ -1329,8 +1594,10 @@ public: ? stm_edge_to_ltm_faces : ltm_edge_to_stm_faces; - compute_intersection_points(tm1_edge_to_tm2_faces, tm1, tm2, vpm1, vpm2, current_node); - compute_intersection_points(tm2_edge_to_tm1_faces, tm2, tm1, vpm2, vpm1, current_node); + compute_intersection_points(tm1_edge_to_tm2_faces, tm1, tm2, vpm1, vpm2, non_manifold_feature_map_1, non_manifold_feature_map_2, current_node); + compute_intersection_points(tm2_edge_to_tm1_faces, tm2, tm1, vpm2, vpm1, non_manifold_feature_map_2, non_manifold_feature_map_1, current_node); + + nodes.check_no_duplicates(); if (!build_polylines){ visitor.finalize(nodes,tm1,tm2,vpm1,vpm2); @@ -1377,7 +1644,7 @@ public: CGAL_assertion(current_node+1==0); //first handle coplanar triangles - compute_intersection_of_coplanar_faces(current_node, tm, tm, vpm, vpm); + compute_intersection_of_coplanar_faces(current_node, tm, tm, vpm, vpm, non_manifold_feature_map_1, non_manifold_feature_map_1); if (!coplanar_faces.empty()) visitor.input_have_coplanar_faces(); @@ -1385,7 +1652,7 @@ public: //compute intersection points of segments and triangles. //build the nodes of the graph and connectivity infos - compute_intersection_points(stm_edge_to_ltm_faces, tm, tm, vpm, vpm, current_node); + compute_intersection_points(stm_edge_to_ltm_faces, tm, tm, vpm, vpm, non_manifold_feature_map_1, non_manifold_feature_map_1, current_node); if (!build_polylines){ visitor.finalize(nodes,tm,tm,vpm,vpm); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h index 99ff0a7815f..ff9c3344330 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h @@ -126,7 +126,13 @@ public: } void all_nodes_created(){} - void finalize() {} + template + void finalize(const Mesh_to_map_node&) {} + + void check_no_duplicates() + { + CGAL_assertion(nodes.size() == std::set(nodes.begin(), nodes.end()).size()); + } }; // end specialization // Intersection_nodes @@ -163,7 +169,6 @@ private: Exact_to_double exact_to_double; Exact_kernel ek; Exact_kernel::Intersect_3 exact_intersection; - std::vector tm1_vertices, tm2_vertices; const bool doing_autorefinement; public: @@ -272,40 +277,46 @@ public: } void all_nodes_created() - { - tm1_vertices.resize(enodes.size(), GT::null_vertex()); - tm2_vertices.resize(enodes.size(), GT::null_vertex()); - } + {} template // VertexPointMap1 or VertexPointMap2 - void call_put(const VPM& vpm, vertex_descriptor vd, std::size_t i, TriangleMesh& tm) + void call_put(const VPM& vpm, vertex_descriptor vd, std::size_t i, TriangleMesh&) { - put(vpm, vd, exact_to_double(enodes[i])); - if (&tm1==&tm) - { - if ( tm1_vertices[i] == GT::null_vertex() ) - { - tm1_vertices[i] = vd; - return; - } - if (doing_autorefinement) - tm2_vertices[i] = vd; - } - else - tm2_vertices[i] = vd; + put(vpm, vd, exact_to_double(enodes[i])); // Note this call is useless and only useful to see something in debug for intermediate results } - void finalize() + template + void finalize(const std::map& mesh_to_node_id_to_vertex) { - for (std::size_t i=0, e=enodes.size(); i!=e; ++i) + if (!doing_autorefinement) { - Point_3 pt = exact_to_double(enodes[i]); - if ( tm1_vertices[i] != GT::null_vertex() ) - put(vpm1, tm1_vertices[i], pt); - if ( tm2_vertices[i] != GT::null_vertex() ) - put(vpm2, tm2_vertices[i], pt); + const Node_id_to_vertex& tm1_vertices = mesh_to_node_id_to_vertex.find(&tm1)->second; + const Node_id_to_vertex& tm2_vertices = mesh_to_node_id_to_vertex.find(&tm2)->second; + for (std::size_t i=0, e=enodes.size(); i!=e; ++i) + { + Point_3 pt = exact_to_double(enodes[i]); + if ( tm1_vertices[i] != GT::null_vertex() ) + put(vpm1, tm1_vertices[i], pt); + if ( tm2_vertices[i] != GT::null_vertex() ) + put(vpm2, tm2_vertices[i], pt); + } + } + else{ + const Node_id_to_vertex& tm1_vertices = mesh_to_node_id_to_vertex.find(&tm1)->second; + for (std::size_t i=0, e=enodes.size(); i!=e; ++i) + { + Point_3 pt = exact_to_double(enodes[i]); + if ( tm1_vertices[i] != GT::null_vertex() ) + put(vpm1, tm1_vertices[i], pt); + } } } + + void check_no_duplicates() + { + CGAL_assertion(enodes.size() == std::set(enodes.begin(), enodes.end()).size()); + } + }; // end specialization // Intersection_nodes @@ -421,9 +432,16 @@ public: } void all_nodes_created(){} - void finalize() {} + + template + void finalize(const std::map&) + {} + void check_no_duplicates() + { + CGAL_assertion(nodes.size() == std::set(nodes.begin(), nodes.end()).size()); + } }; // end specialization diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h index 60877681330..78d96b4c964 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h @@ -1771,6 +1771,11 @@ surface_intersection(const TriangleMesh& tm1, Corefinement::Intersection_of_triangle_meshes functor(tm1, tm2, vpm1, vpm2); + + // Fill non-manifold feature maps if provided + functor.set_non_manifold_feature_map_1(parameters::get_parameter(np1, internal_np::non_manifold_feature_map)); + functor.set_non_manifold_feature_map_2(parameters::get_parameter(np2, internal_np::non_manifold_feature_map)); + return functor(polyline_output, throw_on_self_intersection, true); }