add function to refine a mesh along an isocurve

Installation/CHANGES.md

@@ -52,6 +52,10 @@ Release date: October 2023
  `CGAL::Simplicial_mesh_cell_base_3`
  have been modified to enable passing a geometric traits and a custom cell base class.
 
+### [Polygon Mesh Processing](https://doc.cgal.org/6.0/Manual/packages.html#PkgPolygonMeshProcessing)
+- added the function `CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel()` that refines a polygon mesh
+ along an isocurve.
+
 [Release 5.6](https://github.com/CGAL/cgal/releases/tag/v5.6)
 -----------
 

Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt

@@ -258,6 +258,7 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage.
 - `CGAL::Polygon_mesh_processing::triangle()`
 - `CGAL::Polygon_mesh_processing::region_growing_of_planes_on_faces()`
 - `CGAL::Polygon_mesh_processing::detect_corners_of_regions()`
+- `CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel()`
 
 \cgalCRPSection{I/O Functions}
 - \link PMP_IO_grp `CGAL::Polygon_mesh_processing::IO::read_polygon_mesh()`\endlink

Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt

@@ -73,6 +73,8 @@ if(TARGET CGAL::Eigen3_support)
  target_link_libraries(delaunay_remeshing_example PUBLIC CGAL::Eigen3_support)
  create_single_source_cgal_program("remesh_almost_planar_patches.cpp")
  target_link_libraries(remesh_almost_planar_patches PUBLIC CGAL::Eigen3_support)
+ create_single_source_cgal_program("geodesic_isolevel_refinement.cpp")
+ target_link_libraries(geodesic_isolevel_refinement PUBLIC CGAL::Eigen3_support)
 else()
  message(STATUS "NOTICE: Examples that use Eigen will not be compiled.")
 endif()

Polygon_mesh_processing/examples/Polygon_mesh_processing/geodesic_isolevel_refinement.cpp

@@ -0,0 +1,72 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h>
+
+#include <CGAL/Polygon_mesh_processing/refine_mesh_at_isolevel.h>
+#include <CGAL/Polygon_mesh_processing/connected_components.h>
+
+#include <fstream>
+#include <iostream>
+
+typedef CGAL::Simple_cartesian<double> Kernel;
+typedef Kernel::Point_3 Point_3;
+typedef CGAL::Surface_mesh<Point_3> Triangle_mesh;
+
+typedef boost::graph_traits<Triangle_mesh>::vertex_descriptor vertex_descriptor;
+typedef boost::graph_traits<Triangle_mesh>::edge_descriptor edge_descriptor;
+
+typedef Triangle_mesh::Property_map<vertex_descriptor,double> Vertex_distance_map;
+typedef CGAL::Heat_method_3::Surface_mesh_geodesic_distances_3<Triangle_mesh> Heat_method;
+
+int main(int argc, char* argv[])
+{
+ const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/elephant.off");
+
+ Triangle_mesh tm;
+ if(!CGAL::IO::read_polygon_mesh(filename, tm) ||
+ CGAL::is_empty(tm) || !CGAL::is_triangle_mesh(tm))
+ {
+ std::cerr << "Invalid input file." << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ // default isovalues for cutting the mesh
+ std::vector<double> isovalues = {0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 10};
+
+ if (argc>2)
+ {
+ isovalues.clear();
+ for (int i=2; i<argc; ++i)
+ isovalues.push_back(atof(argv[i]));
+ }
+
+ //property map for the distance values to the source set
+ Vertex_distance_map vertex_distance = tm.add_property_map<vertex_descriptor, double>("v:distance", 0).first;
+
+ Heat_method hm(tm);
+
+ //use heat method to compute approximated geodesic distances to the source vertex `s`
+ vertex_descriptor s = *(vertices(tm).first);
+ hm.add_source(s);
+ hm.estimate_geodesic_distances(vertex_distance);
+
+ // property map to flag new cut edge added in the mesh
+ auto ecm = tm.add_property_map<edge_descriptor, bool>("e:is_constrained", 0).first;
+
+ // refine the mesh along isovalues
+ for (double isovalue : isovalues)
+ CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel(tm, vertex_distance, isovalue, CGAL::parameters::edge_is_constrained_map(ecm));
+
+ // split the mesh in connected components bounded by the isocurves
+ std::vector<Triangle_mesh> edges_split;
+ CGAL::Polygon_mesh_processing::split_connected_components(tm, edges_split, CGAL::parameters::edge_is_constrained_map(ecm));
+
+ assert(argc!=1 || edges_split.size() == 22);
+
+ // export each submesh in a file
+ for(std::size_t i=0; i<edges_split.size(); ++i)
+ std::ofstream("out_"+std::to_string(i)+".off") << std::setprecision(17) << edges_split[i];
+
+ return 0;
+}

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine_mesh_at_isolevel.h

@@ -0,0 +1,174 @@
+// Copyright (c) 2021-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_REFINE_MESH_AT_ISOLEVEL_H
+#define CGAL_POLYGON_MESH_PROCESSING_REFINE_MESH_AT_ISOLEVEL_H
+
+#include <CGAL/license/Polygon_mesh_processing/miscellaneous.h>
+
+#include <CGAL/boost/graph/Euler_operations.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+#include <unordered_set>
+#include <unordered_map>
+#include <unordered_map>
+
+namespace CGAL {
+namespace Polygon_mesh_processing {
+
+/*!
+ * \ingroup PkgPolygonMeshProcessingRef
+ *
+ * refines `pm` by adding new vertices on edges having their incident vertices associated with
+ * values respectively larger and smaller than `isovalue` in `value_map`.
+ * The placement of new vertices on edges will be done by linear interpolation
+ * using the aforementioned values.
+ * New vertices will be associated `isovalue` in `value_map` when created.
+ * Additionally, new edges will be added by connecting new vertices created sharing
+ * a common incident face. Note that in case more than two new vertices are added
+ * on a face boundary, no edges will be created in that face.
+ *
+ * @tparam PolygonMesh a model of the concepts `EdgeListGraph` and `FaceListGraph`
+ * @tparam ValueMap a model of the concept `ReadWritePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
+ * as key type and with its value type being the type of the coordinates of points associated with vertices
+ * in the vertex map provided to the `vertex_point_map()` named parameter.
+ * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" for `pm`
+ *
+ * @param pm the polygon mesh to be refined.
+ * @param value_map the property map containing a value at each vertex for a given function defined over the mesh.
+ * @param isovalue the value used to refine
+ * @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
+ *
+ * \cgalNamedParamsBegin
+ * \cgalParamNBegin{edge_is_constrained_map}
+ * \cgalParamDescription{an output property map associating `true` to all edges connecting vertices on the isolevel,
+ * and `false` for all other edges.}
+ * \cgalParamType{a class model of `WritablePropertyMap` with `boost::graph_traits<PolygonMesh>::%edge_descriptor`
+ * as key type and `bool` as value type}
+ * \cgalParamDefault{No marks on edges will be put}
+ * \cgalParamNEnd
+ *
+ * \cgalParamNBegin{vertex_point_map}
+ * \cgalParamDescription{a property map associating points to the vertices of `pm`}
+ * \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
+ * as key type and `%Point_3` as value type}
+ * \cgalParamDefault{`boost::get(CGAL::vertex_point, pm)`}
+ * \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
+ * must be available in `PolygonMesh`.}
+ * \cgalParamNEnd
+ * \cgalNamedParamsEnd
+ *
+ */
+template <class PolygonMesh, class ValueMap, class NamedParameters = parameters::Default_named_parameters>
+void refine_mesh_at_isolevel(PolygonMesh& pm,
+ ValueMap value_map,
+ typename boost::property_traits<ValueMap>::value_type isovalue,
+ const NamedParameters& np = parameters::default_values())
+{
+ typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
+ typedef typename boost::graph_traits<PolygonMesh>::edge_descriptor edge_descriptor;
+ typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
+ typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
+ typedef typename boost::property_traits<ValueMap>::value_type FT;
+
+ using parameters::choose_parameter;
+ using parameters::get_parameter;
+ using parameters::is_default_parameter;
+
+ typedef Static_boolean_property_map<edge_descriptor, false> Default_ECM;
+ typedef typename internal_np::Lookup_named_param_def<internal_np::edge_is_constrained_t,
+ NamedParameters,
+ Default_ECM>::type ECM;
+ typedef typename GetVertexPointMap < PolygonMesh, NamedParameters>::type VPM;
+
+ VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+ get_property_map(vertex_point, pm));
+
+ ECM ecm = choose_parameter(get_parameter(np, internal_np::edge_is_constrained), Default_ECM());
+
+ std::unordered_map<face_descriptor, std::vector<halfedge_descriptor> > faces_to_split;
+ std::vector<edge_descriptor> to_split;
+ std::unordered_set<vertex_descriptor> vertices_on_isoline;
+ for (edge_descriptor e : edges(pm))
+ {
+ vertex_descriptor src = source(e, pm), tgt = target(e, pm);
+
+ if (get(value_map, src)==isovalue)
+ {
+ vertices_on_isoline.insert(src);
+ if (get(value_map, tgt)==isovalue)
+ {
+ put(ecm, e, true); // special case for faces entirely on an isovalue
+ continue;
+ }
+ continue;
+ }
+ if (get(value_map, tgt)==isovalue)
+ {
+ vertices_on_isoline.insert(tgt);
+ continue;
+ }
+ if ( (get(value_map, tgt) < isovalue) != (get(value_map, src) < isovalue) )
+ {
+ to_split.push_back(e);
+ }
+ }
+
+ for (edge_descriptor e : to_split)
+ {
+ vertex_descriptor src = source(e, pm), tgt = target(e, pm);
+ FT ds = get(value_map, src);
+ FT dt = get(value_map, tgt);
+ FT alpha = (isovalue - dt) / (ds - dt);
+ halfedge_descriptor hnew = CGAL::Euler::split_edge(halfedge(e, pm), pm);
+ put(vpm, target(hnew, pm), barycenter(get(vpm,src), alpha, get(vpm, tgt), 1-alpha));
+ put(value_map, target(hnew, pm) , isovalue);
+ face_descriptor f = face(hnew, pm);
+ if (f!=boost::graph_traits<PolygonMesh>::null_face())
+ faces_to_split[f].push_back(hnew);
+ hnew=pm.prev(opposite(hnew, pm));
+ f = face(hnew, pm);
+ if (f!=boost::graph_traits<PolygonMesh>::null_face())
+ faces_to_split[f].push_back(hnew);
+ }
+
+ for (vertex_descriptor vh : vertices_on_isoline)
+ {
+ for (halfedge_descriptor h : halfedges_around_target(vh, pm))
+ {
+ face_descriptor f = face(h, pm);
+ if (f!=boost::graph_traits<PolygonMesh>::null_face())
+ faces_to_split[f].push_back(h);
+ }
+ }
+
+ for (const auto& p : faces_to_split)
+ {
+ if(p.second.size()!=2) continue;
+
+ std::pair<edge_descriptor ,bool> res = edge(target(p.second[0],pm),
+ target(p.second[1],pm), pm);
+ if (res.second)
+ {
+ // no split as the edge already exists (the two vertices are on the isolevel)
+ put(ecm, res.first, true);
+ continue;
+ }
+
+ halfedge_descriptor hnew = CGAL::Euler::split_face(p.second[0], p.second[1], pm);
+ put(ecm, edge(hnew, pm), true);
+ }
+}
+
+} } // end of CGAL::Polygon_mesh_processing
+
+
+#endif

Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt

@@ -65,6 +65,7 @@ create_single_source_cgal_program("triangulate_hole_with_cdt_2_test.cpp")
 create_single_source_cgal_program("test_pmp_polyhedral_envelope.cpp")
 create_single_source_cgal_program("test_pmp_np_function.cpp")
 create_single_source_cgal_program("test_degenerate_pmp_clip_split_corefine.cpp")
+create_single_source_cgal_program("test_isolevel_refinement.cpp")
 # create_single_source_cgal_program("test_pmp_repair_self_intersections.cpp")
 
 find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater)

Polygon_mesh_processing/test/Polygon_mesh_processing/test_isolevel_refinement.cpp

@@ -0,0 +1,68 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/refine_mesh_at_isolevel.h>
+#include <CGAL/Polygon_mesh_processing/connected_components.h>
+
+#include <fstream>
+#include <iostream>
+
+typedef CGAL::Simple_cartesian<double> Kernel;
+typedef Kernel::Point_3 Point_3;
+typedef CGAL::Surface_mesh<Point_3> Triangle_mesh;
+
+typedef boost::graph_traits<Triangle_mesh>::vertex_descriptor vertex_descriptor;
+typedef boost::graph_traits<Triangle_mesh>::edge_descriptor edge_descriptor;
+
+typedef Triangle_mesh::Property_map<vertex_descriptor,double> Vertex_distance_map;
+
+int main()
+{
+ const std::string filename = CGAL::data_file_path("meshes/elephant.off");
+
+ Triangle_mesh tm;
+ if(!CGAL::IO::read_polygon_mesh(filename, tm) ||
+ CGAL::is_empty(tm) || !CGAL::is_triangle_mesh(tm))
+ {
+ std::cerr << "Invalid input file." << std::endl;
+ return EXIT_FAILURE;
+ }
+ //property map for the distance values to the source set
+ Vertex_distance_map vertex_distance = tm.add_property_map<vertex_descriptor, double>("v:distance", 5).first;
+ std::vector<int> zero_vids={
+ /*a closed polyline*/ 2144,145,2690,1752,339,215,1395,338,77,2145,2052,2054,343,1936,22,1751,214,1499,142,358,2694,1750,301,65,59,2650,2060,205,2651,2061,2490,1939,898,13,298,
+ /*two adjacent triangles*/ 532, 185, 534, 2735,
+ /*another patch with missing crossed edges*/134,73,1883,2533,72,532,185,131,534
+ };
+
+ std::vector<int> minus_vids = {132, 364};
+
+ for (int i : zero_vids)
+ put(vertex_distance, Triangle_mesh::Vertex_index(i), 0);
+ for (int i : minus_vids)
+ put(vertex_distance, Triangle_mesh::Vertex_index(i), -5);
+
+ // property map to flag new cut edge added in the mesh
+ auto ecm = tm.add_property_map<edge_descriptor, bool>("e:is_constrained", 0).first;
+
+ CGAL::Polygon_mesh_processing::refine_mesh_at_isolevel(tm, vertex_distance, 0, CGAL::parameters::edge_is_constrained_map(ecm));
+
+ std::ofstream debug("edges.polylines.txt");
+ for (Triangle_mesh::Edge_index e : edges(tm))
+ if (get(ecm, e))
+ debug << "2 " << tm.point(source(e, tm)) << " " << tm.point(target(e, tm)) << "\n";
+ debug.close();
+
+ // split the mesh in connected components bounded by the isocurves
+ std::vector<Triangle_mesh> edges_split;
+ CGAL::Polygon_mesh_processing::split_connected_components(tm, edges_split, CGAL::parameters::edge_is_constrained_map(ecm));
+
+ // export each submesh in a file
+ for(std::size_t i=0; i<edges_split.size(); ++i)
+ std::ofstream("out_"+std::to_string(i)+".off") << std::setprecision(17) << edges_split[i];
+
+ assert(edges_split.size() == 5);
+
+
+ return 0;
+}