From a1ff847b6a4ec05ce40e3e24b7d38d8a87fd868e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Mon, 2 Jan 2023 08:17:16 +0100 Subject: [PATCH] add "brute-force" version of autorefine() --- .../Polygon_mesh_processing/CMakeLists.txt | 1 + .../triangle_mesh_autorefinement.cpp | 29 + .../Polygon_mesh_processing/autorefinement.h | 507 ++++++++++++++++++ .../internal/parameters_interface.h | 1 + 4 files changed, 538 insertions(+) create mode 100644 Polygon_mesh_processing/examples/Polygon_mesh_processing/triangle_mesh_autorefinement.cpp create mode 100644 Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index 192e7d29b74..e188e01e3eb 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -50,6 +50,7 @@ create_single_source_cgal_program("match_faces.cpp") create_single_source_cgal_program("cc_compatible_orientations.cpp") create_single_source_cgal_program("hausdorff_distance_remeshing_example.cpp") create_single_source_cgal_program("hausdorff_bounded_error_distance_example.cpp") +create_single_source_cgal_program("triangle_mesh_autorefinement.cpp") find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater) include(CGAL_Eigen3_support) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangle_mesh_autorefinement.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangle_mesh_autorefinement.cpp new file mode 100644 index 00000000000..65db52fef2c --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangle_mesh_autorefinement.cpp @@ -0,0 +1,29 @@ +#include +#include + +#include +#include + +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel::Point_3 Point; + +typedef CGAL::Surface_mesh Mesh; + +namespace PMP = CGAL::Polygon_mesh_processing; + +int main(int argc, char** argv) +{ + Mesh mesh; + + const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off") + : std::string(argv[1]); + CGAL::IO::read_polygon_mesh(filename, mesh); + + PMP::autorefine(mesh); + + CGAL::IO::write_polygon_mesh("autorefined.off", mesh, CGAL::parameters::stream_precision(17)); + + return 0; +} diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h new file mode 100644 index 00000000000..8c34b4041d8 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -0,0 +1,507 @@ +// Copyright (c) 2023 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_AUTOREFINEMENT_H +#define CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_H + +#include + +#include +#include +#include +#include + +// output +#include +#include + +#ifndef NDEBUG +// debug +#include +#endif + +#include + +namespace CGAL { +namespace Polygon_mesh_processing { + +#ifndef DOXYGEN_RUNNING +namespace autorefine_impl { + +template +void generate_subtriangles(const typename EK::Triangle_3& t, + const std::vector& segments, + const std::vector& points, + std::vector& new_triangles) +{ + typedef CGAL::Projection_traits_3 P_traits; + // typedef CGAL::Exact_intersections_tag Itag; + typedef CGAL::No_constraint_intersection_requiring_constructions_tag Itag; + typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; + + + P_traits cdt_traits(normal(t[0], t[1], t[2])); + CDT cdt(cdt_traits); + + cdt.insert_outside_affine_hull(t[0]); + cdt.insert_outside_affine_hull(t[1]); + typename CDT::Vertex_handle v = cdt.tds().insert_dim_up(cdt.infinite_vertex(), false); + v->set_point(t[2]); + + for (const typename EK::Segment_3& s : segments) + cdt.insert_constraint(s[0], s[1]); + + for (const typename EK::Point_3& p : points) + cdt.insert(p); + + for (typename CDT::Face_handle fh : cdt.finite_face_handles()) + { + new_triangles.emplace_back(fh->vertex(0)->point(), + fh->vertex(cdt.ccw(0))->point(), + fh->vertex(cdt.cw(0))->point()); + } +} + +template +struct Intersection_visitor +{ + std::vector& all_segments_1; + std::vector& all_segments_2; + std::vector& all_points_1; + std::vector& all_points_2; + + Intersection_visitor(std::vector& all_segments_1, + std::vector& all_segments_2, + std::vector& all_points_1, + std::vector& all_points_2) + : all_segments_1(all_segments_1) + , all_segments_2(all_segments_2) + , all_points_1(all_points_1) + , all_points_2(all_points_2) + {} + + typedef void result_type; + void operator()(const typename EK::Point_3& p) + { + all_points_1.push_back(p); + all_points_2.push_back(p); + } + + void operator()(const typename EK::Segment_3& s) + { + all_segments_1.push_back(s); + all_segments_2.push_back(s); + } + + void operator()(const typename EK::Triangle_3& t) + { + for (std::size_t i=1; i<3; ++i) + { + typename EK::Segment_3 s(t[i-1], t[i]); + all_segments_1.push_back(s); + all_segments_2.push_back(s); + } + } + + void operator()(const std::vector& poly) + { + std::size_t nbp = poly.size(); + for (std::size_t i=1; i +bool is_output_valid(TriangleMesh& tm , VPM vpm, TID_Map tid_map, const std::vector< std::vector>>& triangles) +{ + typedef typename Kernel_traits::value_type>::type IK; + typedef boost::graph_traits Graph_traits; + typedef typename Graph_traits::face_descriptor face_descriptor; + typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor; + + std::vector soup_points; + std::vector > soup_triangles; + Cartesian_converter to_exact; + std::map point_id_map; + + auto get_point_id = [&](const typename EK::Point_3& pt) + { + auto insert_res = point_id_map.insert(std::make_pair(pt, soup_points.size())); + if (insert_res.second) + soup_points.push_back(pt); + return insert_res.first->second; + }; + + for (face_descriptor f : faces(tm)) + { + int tid = get(tid_map, f); + if (tid == -1) + { + halfedge_descriptor h = halfedge(f, tm); + soup_triangles.emplace_back( + CGAL::make_array(get_point_id(to_exact(get(vpm,source(h, tm)))), + get_point_id(to_exact(get(vpm,target(h, tm)))), + get_point_id(to_exact(get(vpm,target(next(h, tm), tm))))) + ); + } + else + { + for (const typename EK::Triangle_3& t : triangles[tid]) + { + soup_triangles.emplace_back(CGAL::make_array(get_point_id(t[0]), get_point_id(t[1]), get_point_id(t[2]))); + } + } + } + + typedef Surface_mesh Exact_mesh; + Exact_mesh etm; + orient_polygon_soup(soup_points, soup_triangles); + polygon_soup_to_polygon_mesh(soup_points, soup_triangles, etm); + typename Exact_mesh::Property_map is_border_map = + etm.template add_property_map("v:is_border", false).first; + for(typename Exact_mesh::Halfedge_index h : etm.halfedges()) + { + if (CGAL::is_border(h, etm)) + is_border_map[target(h, etm)] = true; + } + + //TODO: double check me + auto skip_faces = [&](const std::pair& p) + { + typename Exact_mesh::Halfedge_index h1 = etm.halfedge(p.first), h2=etm.halfedge(p.second); + + boost::container::small_vector bv1; + if (is_border_map[source(h1, etm)]) bv1.push_back(prev(h1, etm)); + if (is_border_map[target(h1, etm)]) bv1.push_back(h1); + if (is_border_map[target(next(h1, etm), etm)]) bv1.push_back(next(h1, etm)); + if (bv1.empty()) return false; + + boost::container::small_vector bv2; + if (is_border_map[source(h2, etm)]) bv2.push_back(prev(h2, etm)); + if (is_border_map[target(h2, etm)]) bv2.push_back(h2); + if (is_border_map[target(next(h2, etm), etm)]) bv2.push_back(next(h2, etm)); + if (bv2.empty()) return false; + + //collect identical border vertices + boost::container::small_vector, 3> common; + for(typename Exact_mesh::Halfedge_index h1 : bv1) + for(typename Exact_mesh::Halfedge_index h2 : bv2) + if (etm.point(target(h1, etm))==etm.point(target(h2,etm))) + common.push_back(std::make_pair(h1,h2)); + + if (common.empty()) return false; + + switch (common.size()) + { + case 1: + { + // geometric check if the opposite segments intersect the triangles + const typename EK::Triangle_3 t1(etm.point(source(h1,etm)), + etm.point(target(h1,etm)), + etm.point(target(next(h1,etm),etm))); + const typename EK::Triangle_3 t2(etm.point(source(h2,etm)), + etm.point(target(h2,etm)), + etm.point(target(next(h2,etm),etm))); + + const typename EK::Segment_3 s1(etm.point(source(common[0].first,etm)), etm.point(target(next(common[0].first,etm),etm))); + const typename EK::Segment_3 s2(etm.point(source(common[0].second,etm)), etm.point(target(next(common[0].second,etm),etm))); + + if(do_intersect(t1, s2) || do_intersect(t2, s1)) + return false; + return true; + } + case 2: + { + // shared edge + h1 = next(common[0].first, etm) == common[1].first ? common[1].first : common[0].first; + h2 = next(common[0].second, etm) == common[1].second ? common[1].second : common[0].second; + + if ( is_border(etm.opposite(h1), etm) && + is_border(etm.opposite(h2), etm) ) + { + if( CGAL::coplanar(etm.point(source(h1,etm)), + etm.point(target(h1,etm)), + etm.point(target(etm.next(h1),etm)), + etm.point(source(h1,etm))) && + CGAL::coplanar_orientation(etm.point(source(h1,etm)), + etm.point(target(h1,etm)), + etm.point(target(etm.next(h1),etm)), + etm.point(source(h1,etm))) == CGAL::POSITIVE) + { + return false; + } + return true; + } + else + { + // TODO: 2 identical border vertices, no common edge on the boundary. Not sure what to do + return false; + } + } + default: // size == 3 + return true; + } + }; + + std::vector< std::pair > si_faces; + self_intersections(etm, + CGAL::filter_output_iterator(std::back_inserter(si_faces), + skip_faces)); + + return si_faces.empty(); +} + +} // end of autorefine_impl +#endif + +/** + * \ingroup PMP_corefinement_grp + * \link coref_def_subsec autorefines \endlink `tm`. Refines a triangle mesh + * so that no triangles intersects in their interior. + * Self-intersection edges will be marked as constrained. If an edge that was marked as + * constrained is split, its sub-edges will be marked as constrained as well. + * + * @tparam TriangleMesh a model of `HalfedgeListGraph`, `FaceListGraph`, and `MutableFaceGraph` + * @tparam NamedParameters a sequence of \ref namedparameters + * + * @param tm input triangulated surface mesh + * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below + * + * \cgalParamNBegin{geom_traits} + * \cgalParamDescription{an instance of a geometric traits class} + * \cgalParamType{a class model of `PMPSelfIntersectionTraits`} + * \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`} + * \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.} + * \cgalParamNEnd + * + * \cgalNamedParamsBegin + * \cgalParamNBegin{vertex_point_map} + * \cgalParamDescription{a property map associating points to the vertices of `tm`} + * \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits::%vertex_descriptor` + * as key type and `%Point_3` as value type} + * \cgalParamDefault{`boost::get(CGAL::vertex_point, tm)`} + * \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t` + * must be available in `TriangleMesh`.} + * \cgalParamNEnd + * + * \cgalParamNBegin{edge_is_constrained_map} + * \cgalParamDescription{a property map containing the constrained-or-not status of each edge of `tm`} + * \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits::%edge_descriptor` + * as key type and `bool` as value type} + * \cgalParamDefault{a constant property map returning `false` for any edge} + * \cgalParamNEnd + * + * \cgalParamNBegin{face_index_map} + * \cgalParamDescription{a property map associating to each face of `tm` a unique index between `0` and `num_faces(tm) - 1`} + * \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits::%face_descriptor` + * as key type and `std::size_t` as value type} + * \cgalParamDefault{an automatically indexed internal map} + * \cgalParamExtra{If the property map is writable, the indices of the faces of `tm1` and `tm2` + * will be set after the corefinement is done.} + * \cgalParamNEnd + * + * \cgalParamNBegin{visitor} + * \cgalParamDescription{a visitor used to track the creation of new faces} + * \cgalParamType{a class model of `PMPCorefinementVisitor`} + * \cgalParamDefault{`Corefinement::Default_visitor`} + * \cgalParamNEnd + * \cgalNamedParamsEnd + * + */ +template +void +autorefine( TriangleMesh& tm, + const NamedParameters& np = parameters::default_values()) +{ + //TODO: what about degenerate faces? + + using parameters::choose_parameter; + using parameters::get_parameter; + + typedef typename GetGeomTraits::type GT; + GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits)); + + typedef typename GetVertexPointMap::type VPM; + VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(vertex_point, tm)); + + typedef typename internal_np::Lookup_named_param_def < + internal_np::concurrency_tag_t, + NamedParameters, + Sequential_tag + > ::type Concurrency_tag; + + typedef boost::graph_traits Graph_traits; + typedef typename Graph_traits::face_descriptor face_descriptor; + typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename Graph_traits::vertex_descriptor vertex_descriptor; + typedef std::pair Pair_of_faces; + + std::vector si_pairs; + + // collect intersecting pairs of triangles + self_intersections(tm, std::back_inserter(si_pairs), np); + + if (si_pairs.empty()) return; + + // assign an id per triangle involved in an intersection + // + the faces involved in the intersection + typedef CGAL::dynamic_face_property_t Face_property_tag; + typedef typename boost::property_map::type Triangle_id_map; + + Triangle_id_map tid_map = get(Face_property_tag(), tm); + for (face_descriptor f : faces(tm)) + put(tid_map, f, -1); + + std::vector intersected_faces; + int tid=-1; + for (const Pair_of_faces& p : si_pairs) + { + if (get(tid_map, p.first)==-1) + { + put(tid_map, p.first, ++tid); + intersected_faces.push_back(p.first); + } + if (get(tid_map, p.second)==-1) + { + put(tid_map, p.second, ++tid); + intersected_faces.push_back(p.second); + } + } + + // init the vector of triangles used for the autorefinement of triangles + typedef CGAL::Exact_predicates_exact_constructions_kernel EK; + std::vector< std::vector > triangles(tid+1); + Cartesian_converter to_exact; + + for(face_descriptor f : intersected_faces) + { + halfedge_descriptor h = halfedge(f, tm); + triangles[get(tid_map, f)].emplace_back( + to_exact( get(vpm, source(h, tm)) ), + to_exact( get(vpm, target(h, tm)) ), + to_exact( get(vpm, target(next(h, tm), tm)) ) ); + } + + typename EK::Intersect_3 intersection = EK().intersect_3_object(); + for (const Pair_of_faces& p : si_pairs) + { + int i1 = get(tid_map, p.first), + i2 = get(tid_map, p.second); + + + std::size_t nbt_1 = triangles[i1].size(), + nbt_2 = triangles[i2].size(); + + std::vector< std::vector > all_segments_1(nbt_1); + std::vector< std::vector > all_segments_2(nbt_2); + std::vector< std::vector > all_points_1(nbt_1); + std::vector< std::vector > all_points_2(nbt_2); + + std::vector t1_subtriangles, t2_subtriangles; + for (std::size_t it1=0; it1 intersection_visitor(all_segments_1[it1], all_segments_2[it2], + all_points_1[it1], all_points_2[it2]); + + boost::apply_visitor(intersection_visitor, *inter); + } + } + } + + // now refine triangles + std::vector new_triangles; + for(std::size_t it1=0; it1(triangles[i1][it1], all_segments_1[it1], all_points_1[it1], new_triangles); + } + triangles[i1].swap(new_triangles); + new_triangles.clear(); + for(std::size_t it2=0; it2(triangles[i2][it2], all_segments_2[it2], all_points_2[it2], new_triangles); + } + triangles[i2].swap(new_triangles); + } + + CGAL_assertion( autorefine_impl::is_output_valid(tm, vpm, tid_map, triangles) ); + + // brute force output: create a soup, orient and to-mesh + // WARNING: there is no reason when using double that identical exact points are identical in double + std::vector soup_points; + std::vector > soup_triangles; + Cartesian_converter to_input; + std::map point_id_map; + + for (vertex_descriptor v : vertices(tm)) + { + if (point_id_map.insert(std::make_pair(to_exact(get(vpm,v)), soup_points.size())).second) + soup_points.push_back(get(vpm,v)); + } + + auto get_point_id = [&](const typename EK::Point_3& pt) + { + auto insert_res = point_id_map.insert(std::make_pair(pt, soup_points.size())); + if (insert_res.second) + soup_points.push_back(to_input(pt)); + return insert_res.first->second; + }; + + for (face_descriptor f : faces(tm)) + { + int tid = get(tid_map, f); + if (tid == -1) + { + halfedge_descriptor h = halfedge(f, tm); + soup_triangles.emplace_back( + CGAL::make_array(get_point_id(to_exact(get(vpm,source(h, tm)))), + get_point_id(to_exact(get(vpm,target(h, tm)))), + get_point_id(to_exact(get(vpm,target(next(h, tm), tm))))) + ); + } + else + { + for (const typename EK::Triangle_3& t : triangles[tid]) + { + soup_triangles.emplace_back(CGAL::make_array(get_point_id(t[0]), get_point_id(t[1]), get_point_id(t[2]))); + } + } + } + clear(tm); + orient_polygon_soup(soup_points, soup_triangles); + polygon_soup_to_polygon_mesh(soup_points, soup_triangles, tm); +} + +} } // end of CGAL::Polygon_mesh_processing + +#endif // CGAL_POLYGON_MESH_PROCESSING_AUTOREFINEMENT_H diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h index 5c649faa56c..7d758abc4ee 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -45,6 +45,7 @@ CGAL_add_named_parameter(implementation_tag_t, implementation_tag, implementatio CGAL_add_named_parameter(prevent_unselection_t, prevent_unselection, prevent_unselection) CGAL_add_named_parameter(verbose_t, verbose, verbose) +CGAL_add_named_parameter(concurrency_tag_t, concurrency_tag, concurrency_tag) // List of named parameters used for IO CGAL_add_named_parameter(vertex_normal_output_iterator_t, vertex_normal_output_iterator, vertex_normal_output_iterator)