diff --git a/Classification/examples/Classification/gis_tutorial_example.cpp b/Classification/examples/Classification/gis_tutorial_example.cpp index 7e214e481f5..e02d128a307 100644 --- a/Classification/examples/Classification/gis_tutorial_example.cpp +++ b/Classification/examples/Classification/gis_tutorial_example.cpp @@ -670,7 +670,7 @@ int main (int argc, char** argv) } std::size_t nb_vertices - = std::accumulate (polylines.begin(), polylines.end(), 0u, + = std::accumulate (polylines.begin(), polylines.end(), std::size_t(0), [](std::size_t size, const std::vector& poly) -> std::size_t { return size + poly.size(); }); diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index 8a9d67a73bd..efd83ad75ec 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -44,6 +44,11 @@ Release date: June 2022 ### [2D Regularized Boolean Set-Operations](https://doc.cgal.org/5.5/Manual/packages.html#PkgBooleanSetOperations2) - The concept `GeneralPolygonSetTraits_2` now requires the nested type `Construct_polygon_with_holes_2` instead of `Construct_general_polygon_with_holes_2`. +### [Polygon Mesh Processing](https://doc.cgal.org/5.5/Manual/packages.html#PkgPolygonMeshProcessing) + +- Added the function `CGAL::Polygon_mesh_processing::tangential_relaxation()`, which +applies an area-based tangential mesh smoothing to the vertices of a surface triangle mesh. + [Release 5.4](https://github.com/CGAL/cgal/releases/tag/v5.4) ----------- diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt index f9a9a284d3b..014ccb58919 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt @@ -61,6 +61,9 @@ /// \defgroup PMP_IO_grp I/O Functions /// \ingroup PkgPolygonMeshProcessingRef +/// \defgroup PMPDeprecated Deprecated Functions +/// \ingroup PkgPolygonMeshProcessingRef + /*! \addtogroup PkgPolygonMeshProcessingRef @@ -111,7 +114,9 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage. - `CGAL::Polygon_mesh_processing::triangulate_face()` - `CGAL::Polygon_mesh_processing::triangulate_faces()` - `CGAL::Polygon_mesh_processing::extrude_mesh()` -- `CGAL::Polygon_mesh_processing::smooth_mesh()` +- `CGAL::Polygon_mesh_processing::smooth_mesh()` (deprecated) +- `CGAL::Polygon_mesh_processing::angle_and_area_smoothing()` +- `CGAL::Polygon_mesh_processing::tangential_relaxation()` - `CGAL::Polygon_mesh_processing::smooth_shape()` - `CGAL::Polygon_mesh_processing::random_perturbation()` diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt index 90c62b0fd66..19e6222a76b 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt @@ -138,7 +138,8 @@ While mesh smoothing is achieved by improving the quality of triangles based on shape smoothing is designed to be \e intrinsic, depending as little as possible on the discretization and smoothing the shape alone without optimizing the shape of the triangles. -- Mesh smoothing: `CGAL::Polygon_mesh_processing::smooth_mesh()` moves vertices to optimize geometry around each vertex: +- Mesh smoothing by angle and area optimization: +`CGAL::Polygon_mesh_processing::angle_and_area_smoothing()` moves vertices to optimize geometry around each vertex: it can try to equalize the angles between incident edges, or (and) move vertices so that areas of adjacent triangles tend to equalize. Border vertices are considered constrained and do not move at any step of the procedure. No vertices are inserted at any time. Angle and area smoothing algorithms are based on Surazhsky and Gotsman \cgalCite{cgal:sg-hqct-04}. @@ -165,8 +166,15 @@ Mesh smoothing of the closed surface blobby, containing self-intersection Statistics for the various combinations of mesh smoothing. \cgalFigureEnd +- Mesh smoothing by tangential relaxation: `CGAL::Polygon_mesh_processing::tangential_relaxation()` +moves vertices following an area-based Laplacian smoothing scheme, performed +at each vertex in an estimated tangent plane to the surface. +The full algorithm is described in \cgalCite{botsch2010PMP}. +The example \ref Polygon_mesh_processing/tangential_relaxation_example.cpp shows how this +mesh relaxation function can be used. + - Shape smoothing: `CGAL::Polygon_mesh_processing::smooth_shape()` -incrementally moves vertices towards a weighted barycenter of their neighbors along the mean curvature flow. +moves vertices towards a weighted barycenter of their neighbors along the mean curvature flow. The curvature flow algorithm for shape smoothing is based on Desbrun et al. \cgalCite{cgal:dmsb-ifamdcf-99} and on Kazhdan et al. \cgalCite{kazhdan2012can}. The algorithm uses the mean curvature flow to calculate the translation of vertices along the surface normal with a speed equal to the mean curvature of the area that is being smoothed. This means that vertices on sharp corners slide faster. diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt index 1b053016791..1e172cdaa69 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt @@ -32,6 +32,7 @@ \example Polygon_mesh_processing/repair_polygon_soup_example.cpp \example Polygon_mesh_processing/mesh_smoothing_example.cpp \example Polygon_mesh_processing/shape_smoothing_example.cpp +\example Polygon_mesh_processing/tangential_relaxation_example.cpp \example Polygon_mesh_processing/locate_example.cpp \example Polygon_mesh_processing/orientation_pipeline_example.cpp \example Polygon_mesh_processing/triangulate_faces_split_visitor_example.cpp diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index 1a70c77eab4..c83b338a588 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -73,6 +73,7 @@ create_single_source_cgal_program( "mesh_slicer_example.cpp") #create_single_source_cgal_program( "remove_degeneracies_example.cpp") create_single_source_cgal_program("isotropic_remeshing_example.cpp") create_single_source_cgal_program("isotropic_remeshing_of_patch_example.cpp") +create_single_source_cgal_program("tangential_relaxation_example.cpp") create_single_source_cgal_program("surface_mesh_intersection.cpp") create_single_source_cgal_program("corefinement_SM.cpp") create_single_source_cgal_program("corefinement_consecutive_bool_op.cpp") diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_smoothing_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_smoothing_example.cpp index d177a2a3781..ac7a4739470 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_smoothing_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_smoothing_example.cpp @@ -1,7 +1,7 @@ #include #include -#include +#include #include #include @@ -43,9 +43,9 @@ int main(int argc, char** argv) std::cout << "Smoothing mesh... (" << nb_iterations << " iterations)" << std::endl; // Smooth with both angle and area criteria + Delaunay flips - PMP::smooth_mesh(mesh, CGAL::parameters::number_of_iterations(nb_iterations) - .use_safety_constraints(false) // authorize all moves - .edge_is_constrained_map(eif)); + PMP::angle_and_area_smoothing(mesh, CGAL::parameters::number_of_iterations(nb_iterations) + .use_safety_constraints(false) // authorize all moves + .edge_is_constrained_map(eif)); CGAL::IO::write_polygon_mesh("mesh_smoothed.off", mesh, CGAL::parameters::stream_precision(17)); diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/tangential_relaxation_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/tangential_relaxation_example.cpp new file mode 100644 index 00000000000..c915e2ea419 --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/tangential_relaxation_example.cpp @@ -0,0 +1,35 @@ +#include +#include +#include +#include + +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Mesh; + + +namespace PMP = CGAL::Polygon_mesh_processing; + + +int main(int argc, char* argv[]) +{ + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/pig.off"); + + Mesh mesh; + if(!PMP::IO::read_polygon_mesh(filename, mesh) || !CGAL::is_triangle_mesh(mesh)) + { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + unsigned int nb_iter = (argc > 2) ? std::stoi(std::string(argv[2])) : 10; + + std::cout << "Relax..."; + + PMP::tangential_relaxation(mesh, CGAL::parameters::number_of_iterations(nb_iter)); + + std::cout << "done." << std::endl; + + return 0; +} diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/angle_and_area_smoothing.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/angle_and_area_smoothing.h new file mode 100644 index 00000000000..b34fb9331cf --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/angle_and_area_smoothing.h @@ -0,0 +1,361 @@ +// Copyright (c) 2018 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) : Mael Rouxel-Labbé +// Konstantinos Katrioplas (konst.katrioplas@gmail.com) + +#ifndef CGAL_POLYGON_MESH_PROCESSING_ANGLE_AND_AREA_SMOOTHING_H +#define CGAL_POLYGON_MESH_PROCESSING_ANGLE_AND_AREA_SMOOTHING_H + +#include + +#include +#include +#include + +#include +#include + +#include + +namespace CGAL { +namespace Polygon_mesh_processing { + +/*! +* \ingroup PMP_meshing_grp +* +* \short smooths a triangulated region of a polygon mesh. +* +* This function attempts to make the triangle angle and area distributions as uniform as possible +* by moving (non-constrained) vertices. +* +* Angle-based smoothing does not change the combinatorial information of the mesh. Area-based smoothing +* might change the combinatorial information, unless specified otherwise. It is also possible +* to make the smoothing algorithm "safer" by rejecting moves that, when applied, would worsen the +* quality of the mesh, e.g. that would decrease the value of the smallest angle around a vertex or +* create self-intersections. +* +* Optionally, the points are reprojected after each iteration. +* +* @tparam TriangleMesh model of `MutableFaceGraph`. +* @tparam FaceRange range of `boost::graph_traits::%face_descriptor`, + model of `Range`. Its iterator type is `ForwardIterator`. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" +* +* @param tmesh a polygon mesh with triangulated surface patches to be smoothed. +* @param faces the range of triangular faces defining one or several surface patches to be smoothed. +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below +* +* \cgalNamedParamsBegin +* \cgalParamNBegin{number_of_iterations} +* \cgalParamDescription{the number of iterations for the sequence of the smoothing iterations performed} +* \cgalParamType{unsigned int} +* \cgalParamDefault{`1`} +* \cgalParamNEnd +* +* \cgalParamNBegin{use_angle_smoothing} +* \cgalParamDescription{value to indicate whether angle-based smoothing should be used} +* \cgalParamType{Boolean} +* \cgalParamDefault{`true`} +* \cgalParamNEnd +* +* \cgalParamNBegin{use_area_smoothing} +* \cgalParamDescription{value to indicate whether area-based smoothing should be used} +* \cgalParamType{Boolean} +* \cgalParamDefault{`true`} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_point_map} +* \cgalParamDescription{a property map associating points to the vertices of `tmesh`} +* \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` +* as key type and `%Point_3` as value type} +* \cgalParamDefault{`boost::get(CGAL::vertex_point, tmesh)`} +* \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t` +* must be available in `TriangleMesh`.} +* \cgalParamNEnd +* +* \cgalParamNBegin{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class} +* \cgalParamType{a class model of `Kernel`} +* \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 +* +* \cgalParamNBegin{use_safety_constraints} +* \cgalParamDescription{If `true`, vertex moves that would worsen the mesh are ignored.} +* \cgalParamType{Boolean} +* \cgalParamDefault{`false`} +* \cgalParamNEnd +* +* \cgalParamNBegin{use_Delaunay_flips} +* \cgalParamDescription{If `true`, area-based smoothing will be completed by a phase of +* Delaunay-based edge-flips to prevent the creation of elongated triangles.} +* \cgalParamType{Boolean} +* \cgalParamDefault{`true`} +* \cgalParamNEnd +* +* \cgalParamNBegin{do_project} +* \cgalParamDescription{If `true`, points are projected onto the initial surface after each iteration.} +* \cgalParamType{Boolean} +* \cgalParamDefault{`true`} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_is_constrained_map} +* \cgalParamDescription{a property map containing the constrained-or-not status of each vertex of `tmesh`.} +* \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` +* as key type and `bool` as value type. It must be default constructible.} +* \cgalParamDefault{a default property map where no vertex is constrained} +* \cgalParamExtra{A constrained vertex cannot be modified at all during smoothing.} +* \cgalParamNEnd +* +* \cgalParamNBegin{edge_is_constrained_map} +* \cgalParamDescription{a property map containing the constrained-or-not status of each edge of `tmesh`.} +* \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits::%edge_descriptor` +* as key type and `bool` as value type. It must be default constructible.} +* \cgalParamDefault{a default property map where no edge is constrained} +* \cgalParamExtra{A constrained edge cannot be modified at all during smoothing.} +* \cgalParamNEnd +* \cgalNamedParamsEnd +* +* @warning The third party library \link thirdpartyCeres Ceres \endlink is required +* to use area-based smoothing. +* +* @pre `tmesh` does not contain any degenerate faces. +*/ +template +void angle_and_area_smoothing(const FaceRange& faces, + TriangleMesh& tmesh, + const NamedParameters& np = parameters::default_values()) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typedef typename GetGeomTraits::type GeomTraits; + typedef typename GetVertexPointMap::type VertexPointMap; + + // We need a default pmap that is not just 'constant_pmap(false)' because if an edge is constrained, + // its vertices are constrained. + typedef CGAL::dynamic_vertex_property_t Vertex_property_tag; + typedef typename boost::property_map::type Default_VCMap; + typedef typename internal_np::Lookup_named_param_def ::type VCMap; + + typedef typename internal_np::Lookup_named_param_def // default + > ::type ECMap; + + typedef internal::Area_smoother Area_optimizer; + typedef internal::Mesh_smoother Area_smoother; + typedef internal::Delaunay_edge_flipper Delaunay_flipper; + + typedef internal::Angle_smoother Angle_optimizer; + typedef internal::Mesh_smoother Angle_smoother; + + typedef typename GeomTraits::Triangle_3 Triangle; + typedef std::vector Triangle_container; + + typedef CGAL::AABB_triangle_primitive AABB_Primitive; + typedef CGAL::AABB_traits AABB_Traits; + typedef CGAL::AABB_tree Tree; + + if(std::begin(faces) == std::end(faces)) + return; + + using parameters::choose_parameter; + using parameters::get_parameter; + + // named parameters + GeomTraits gt = choose_parameter(get_parameter(np, internal_np::geom_traits)); + VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(CGAL::vertex_point, tmesh)); + + const bool use_angle_smoothing = choose_parameter(get_parameter(np, internal_np::use_angle_smoothing), true); + bool use_area_smoothing = choose_parameter(get_parameter(np, internal_np::use_area_smoothing), true); + +#ifndef CGAL_PMP_USE_CERES_SOLVER + std::cerr << "Area-based smoothing requires the Ceres Library, which is not available." << std::endl; + std::cerr << "No such smoothing will be performed!" << std::endl; + use_area_smoothing = false; +#endif + + if(!use_angle_smoothing && !use_area_smoothing) + std::cerr << "Called PMP::angle_and_area_smoothing() without any smoothing method selected or available" << std::endl; + + unsigned int nb_iterations = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 1); + const bool do_project = choose_parameter(get_parameter(np, internal_np::do_project), true); + const bool use_safety_constraints = choose_parameter(get_parameter(np, internal_np::use_safety_constraints), true); + const bool use_Delaunay_flips = choose_parameter(get_parameter(np, internal_np::use_Delaunay_flips), true); + + VCMap vcmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained), + get(Vertex_property_tag(), tmesh)); + + // If it's the default vcmap, manually set everything to false because the dynamic pmap has no default initialization + if((std::is_same::value)) + { + for(vertex_descriptor v : vertices(tmesh)) + put(vcmap, v, false); + } + + ECMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained), + Static_boolean_property_map()); + + // a constrained edge has constrained extremities + for(face_descriptor f : faces) + { + if(f == boost::graph_traits::null_face()) + continue; + + for(halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(f, tmesh), tmesh)) + { + if(get(ecmap, edge(h, tmesh))) + { + put(vcmap, source(h, tmesh), true); + put(vcmap, target(h, tmesh), true); + } + } + } + + // Construct the AABB tree (if needed for reprojection) + std::vector input_triangles; + + if(do_project) + { + input_triangles.reserve(faces.size()); + + for(face_descriptor f : faces) + { + halfedge_descriptor h = halfedge(f, tmesh); + if(is_border(h, tmesh)) // should not happen, but just in case + continue; + + input_triangles.push_back(gt.construct_triangle_3_object()(get(vpmap, source(h, tmesh)), + get(vpmap, target(h, tmesh)), + get(vpmap, target(next(h, tmesh), tmesh)))); + } + } + + Tree aabb_tree(input_triangles.begin(), input_triangles.end()); + + // Setup the working ranges and check some preconditions + Angle_smoother angle_smoother(tmesh, vpmap, vcmap, gt); + Area_smoother area_smoother(tmesh, vpmap, vcmap, gt); + Delaunay_flipper delaunay_flipper(tmesh, vpmap, ecmap, gt); + + if(use_angle_smoothing) + angle_smoother.init_smoothing(faces); + + if(use_area_smoothing) + area_smoother.init_smoothing(faces); + + for(unsigned int i=0; i +void angle_and_area_smoothing(TriangleMesh& tmesh, const CGAL_NP_CLASS& np = parameters::default_values()) +{ + angle_and_area_smoothing(faces(tmesh), tmesh, np); +} + +template +void angles_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output) +{ + internal::Quality_evaluator evaluator(tmesh, traits); + evaluator.gather_angles(); + evaluator.extract_angles(output); +} + +template +void areas_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output) +{ + internal::Quality_evaluator evaluator(tmesh, traits); + evaluator.measure_areas(); + evaluator.extract_areas(output); +} + +template +void aspect_ratio_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output) +{ + internal::Quality_evaluator evaluator(tmesh, traits); + evaluator.calc_aspect_ratios(); + evaluator.extract_aspect_ratios(output); +} +///\endcond SKIP_IN_MANUAL + +} // namespace Polygon_mesh_processing +} // namespace CGAL + +#endif // CGAL_POLYGON_MESH_PROCESSING_ANGLE_AND_AREA_SMOOTHING_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index 8c5cb13c0ae..12d0a273651 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h @@ -22,6 +22,7 @@ #include #include #include +#include #include #include @@ -43,6 +44,7 @@ #include #include #include +#include #include #include @@ -66,6 +68,7 @@ #ifdef CGAL_PMP_REMESHING_VERBOSE_PROGRESS #define CGAL_PMP_REMESHING_VERBOSE +#define CGAL_PMP_TANGENTIAL_RELAXATION_VERBOSE #endif namespace CGAL { @@ -980,107 +983,50 @@ namespace internal { // "applies an iterative smoothing filter to the mesh. // The vertex movement has to be constrained to the vertex tangent plane [...] // smoothing algorithm with uniform Laplacian weights" - void tangential_relaxation(const bool relax_constraints/*1d smoothing*/ - , const unsigned int nb_iterations) + void tangential_relaxation_impl(const bool relax_constraints/*1d smoothing*/ + , const unsigned int nb_iterations) { #ifdef CGAL_PMP_REMESHING_VERBOSE std::cout << "Tangential relaxation (" << nb_iterations << " iter.)..."; std::cout << std::endl; #endif - for (unsigned int nit = 0; nit < nb_iterations; ++nit) + + // property map of constrained edges for relaxation + auto edge_constraint = [&](const edge_descriptor e) { -#ifdef CGAL_PMP_REMESHING_VERBOSE_PROGRESS - std::cout << "\r\t(iteration " << (nit + 1) << " / "; - std::cout << nb_iterations << ") "; - std::cout.flush(); -#endif - typedef std::tuple VNP; - std::vector< VNP > barycenters; - // at each vertex, compute vertex normal - // at each vertex, compute barycenter of neighbors - for(vertex_descriptor v : vertices(mesh_)) + return this->is_constrained(e); + }; + auto constrained_edges_pmap + = boost::make_function_property_map(edge_constraint); + + // property map of constrained vertices for relaxation + auto vertex_constraint = [&](const vertex_descriptor v) + { + for (halfedge_descriptor h : halfedges_around_target(v, mesh_)) { - if (is_constrained(v) || is_isolated(v)) - continue; - - else if (is_on_patch(v)) - { - Vector_3 vn = compute_vertex_normal(v, mesh_, parameters::vertex_point_map(vpmap_).geom_traits(gt_)); - Vector_3 move = CGAL::NULL_VECTOR; - unsigned int star_size = 0; - for(halfedge_descriptor h : halfedges_around_target(v, mesh_)) - { - move = move + Vector_3(get(vpmap_, v), get(vpmap_, source(h, mesh_))); - ++star_size; - } - CGAL_assertion(star_size > 0); //isolated vertices have already been discarded - move = (1. / (double)star_size) * move; - - barycenters.push_back( VNP(v, vn, get(vpmap_, v) + move) ); - } - else if (relax_constraints - && !protect_constraints_ - && is_on_patch_border(v) - && !is_corner(v)) - { - Vector_3 vn(NULL_VECTOR); - - std::vector border_halfedges; - for(halfedge_descriptor h : halfedges_around_target(v, mesh_)) - { - if (is_on_patch_border(h) || is_on_patch_border(opposite(h, mesh_))) - border_halfedges.push_back(h); - } - if (border_halfedges.size() == 2)//others are corner cases - { - vertex_descriptor ph0 = source(border_halfedges[0], mesh_); - vertex_descriptor ph1 = source(border_halfedges[1], mesh_); - double dot = to_double(Vector_3(get(vpmap_, v), get(vpmap_, ph0)) - * Vector_3(get(vpmap_, v), get(vpmap_, ph1))); - //check squared cosine is < 0.25 (~120 degrees) - if (0.25 < dot / (sqlength(border_halfedges[0]) * sqlength(border_halfedges[0]))) - barycenters.push_back( VNP(v, vn, CGAL::midpoint(midpoint(border_halfedges[0]), - midpoint(border_halfedges[1]))) ); - } - } + Halfedge_status s = status(h); + if ( s == PATCH + || s == PATCH_BORDER + || status(opposite(h, mesh_)) == PATCH_BORDER) + return false; } + return true; + }; + auto constrained_vertices_pmap + = boost::make_function_property_map(vertex_constraint); - // compute moves - typedef std::pair VP_pair; - std::vector< std::pair > new_locations; - new_locations.reserve(barycenters.size()); - for(const VNP& vnp : barycenters) - { - vertex_descriptor v = std::get<0>(vnp); - Point pv = get(vpmap_, v); - const Vector_3& nv = std::get<1>(vnp); - const Point& qv = std::get<2>(vnp); //barycenter at v + tangential_relaxation( + vertices(mesh_), + mesh_, + CGAL::parameters::number_of_iterations(nb_iterations) + .vertex_point_map(vpmap_) + .geom_traits(gt_) + .edge_is_constrained_map(constrained_edges_pmap) + .vertex_is_constrained_map(constrained_vertices_pmap) + .relax_constraints(relax_constraints) + ); - new_locations.push_back( std::make_pair(v, qv + (nv * Vector_3(qv, pv)) * nv) ); - } - - // perform moves - for(const VP_pair& vp : new_locations) - { - const Point initial_pos = get(vpmap_, vp.first); - const Vector_3 move(initial_pos, vp.second); - - put(vpmap_, vp.first, vp.second); - - //check that no inversion happened - double frac = 1.; - while (frac > 0.03 //5 attempts maximum - && !check_normals(vp.first)) //if a face has been inverted - { - frac = 0.5 * frac; - put(vpmap_, vp.first, initial_pos + frac * move);//shorten the move by 2 - } - if (frac <= 0.02) - put(vpmap_, vp.first, initial_pos);//cancel move - } - - CGAL_assertion(!input_mesh_is_valid_ || is_valid_polygon_mesh(mesh_)); - }//end for loop (nit == nb_iterations) + CGAL_assertion(!input_mesh_is_valid_ || is_valid_polygon_mesh(mesh_)); #ifdef CGAL_PMP_REMESHING_DEBUG debug_self_intersections(); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h index 9e4a5912c82..904bb3f0659 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h @@ -160,7 +160,7 @@ namespace Polygon_mesh_processing { * * \cgalParamNBegin{relax_constraints} * \cgalParamDescription{If `true`, the end vertices of the edges set as constrained -* in `edge_is_constrained_map` and boundary edges move along the} +* in `edge_is_constrained_map` and boundary edges move along the * constrained polylines they belong to.} * \cgalParamType{Boolean} * \cgalParamDefault{`false`} @@ -322,7 +322,7 @@ void isotropic_remeshing(const FaceRange& faces } if(do_flip) remesher.flip_edges_for_valence_and_shape(); - remesher.tangential_relaxation(smoothing_1d, nb_laplacian); + remesher.tangential_relaxation_impl(smoothing_1d, nb_laplacian); if ( choose_parameter(get_parameter(np, internal_np::do_project), true) ) remesher.project_to_surface(get_parameter(np, internal_np::projection_functor)); #ifdef CGAL_PMP_REMESHING_VERBOSE diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h index 87c3bc18fd3..b69aab6dd91 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h @@ -22,7 +22,7 @@ #include #include #include -#include +#include #include #ifndef CGAL_PMP_REMOVE_SELF_INTERSECTION_NO_POLYHEDRAL_ENVELOPE_CHECK #include @@ -546,9 +546,11 @@ bool remove_self_intersections_with_smoothing(std::set -#include -#include -#include +#define CGAL_DEPRECATED_HEADER "" +#define CGAL_REPLACEMENT_HEADER "" +#include -#include -#include +#include -#include +#ifndef CGAL_NO_DEPRECATED_CODE namespace CGAL { namespace Polygon_mesh_processing { /*! -* \ingroup PMP_meshing_grp +* \ingroup PMPDeprecated * -* \brief smooths a triangulated region of a polygon mesh. -* -* This function attempts to make the triangle angle and area distributions as uniform as possible -* by moving (non-constrained) vertices. -* -* Angle-based smoothing does not change the combinatorial information of the mesh. Area-based smoothing -* might change the combinatorial information, unless specified otherwise. It is also possible -* to make the smoothing algorithm "safer" by rejecting moves that, when applied, would worsen the -* quality of the mesh, e.g. that would decrease the value of the smallest angle around a vertex or -* create self-intersections. -* -* Optionally, the points are reprojected after each iteration. -* -* @tparam TriangleMesh model of `MutableFaceGraph`. -* @tparam FaceRange range of `boost::graph_traits::%face_descriptor`, - model of `Range`. Its iterator type is `ForwardIterator`. -* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" -* -* @param tmesh a polygon mesh with triangulated surface patches to be smoothed. -* @param faces the range of triangular faces defining one or several surface patches to be smoothed. -* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below -* -* \cgalNamedParamsBegin -* \cgalParamNBegin{number_of_iterations} -* \cgalParamDescription{the number of iterations for the sequence of the smoothing iterations performed} -* \cgalParamType{unsigned int} -* \cgalParamDefault{`1`} -* \cgalParamNEnd -* -* \cgalParamNBegin{use_angle_smoothing} -* \cgalParamDescription{value to indicate whether angle-based smoothing should be used} -* \cgalParamType{Boolean} -* \cgalParamDefault{`true`} -* \cgalParamNEnd -* -* \cgalParamNBegin{use_area_smoothing} -* \cgalParamDescription{value to indicate whether area-based smoothing should be used} -* \cgalParamType{Boolean} -* \cgalParamDefault{`true`} -* \cgalParamNEnd -* -* \cgalParamNBegin{vertex_point_map} -* \cgalParamDescription{a property map associating points to the vertices of `tmesh`} -* \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` -* as key type and `%Point_3` as value type} -* \cgalParamDefault{`boost::get(CGAL::vertex_point, tmesh)`} -* \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t` -* must be available in `TriangleMesh`.} -* \cgalParamNEnd -* -* \cgalParamNBegin{geom_traits} -* \cgalParamDescription{an instance of a geometric traits class} -* \cgalParamType{a class model of `Kernel`} -* \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 -* -* \cgalParamNBegin{use_safety_constraints} -* \cgalParamDescription{If `true`, vertex moves that would worsen the mesh are ignored.} -* \cgalParamType{Boolean} -* \cgalParamDefault{`false`} -* \cgalParamNEnd -* -* \cgalParamNBegin{use_Delaunay_flips} -* \cgalParamDescription{If `true`, area-based smoothing will be completed by a phase of -* Delaunay-based edge-flips to prevent the creation of elongated triangles.} -* \cgalParamType{Boolean} -* \cgalParamDefault{`true`} -* \cgalParamNEnd -* -* \cgalParamNBegin{do_project} -* \cgalParamDescription{If `true`, points are projected onto the initial surface after each iteration.} -* \cgalParamType{Boolean} -* \cgalParamDefault{`true`} -* \cgalParamNEnd -* -* \cgalParamNBegin{vertex_is_constrained_map} -* \cgalParamDescription{a property map containing the constrained-or-not status of each vertex of `tmesh`.} -* \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` -* as key type and `bool` as value type. It must be default constructible.} -* \cgalParamDefault{a default property map where no vertex is constrained} -* \cgalParamExtra{A constrained vertex cannot be modified at all during smoothing.} -* \cgalParamNEnd -* -* \cgalParamNBegin{edge_is_constrained_map} -* \cgalParamDescription{a property map containing the constrained-or-not status of each edge of `tmesh`.} -* \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits::%edge_descriptor` -* as key type and `bool` as value type. It must be default constructible.} -* \cgalParamDefault{a default property map where no edge is constrained} -* \cgalParamExtra{A constrained edge cannot be modified at all during smoothing.} -* \cgalParamNEnd -* \cgalNamedParamsEnd -* -* @warning The third party library \link thirdpartyCeres Ceres \endlink is required -* to use area-based smoothing. -* -* @pre `tmesh` does not contain any degenerate faces -* -* @see `smooth_shape()` +* \deprecated This function is deprecated since \cgal 5.5, +* `CGAL::angle_and_area_smoothing()` should be used instead. */ template -void smooth_mesh(const FaceRange& faces, - TriangleMesh& tmesh, - const NamedParameters& np = parameters::default_values()) +CGAL_DEPRECATED void smooth_mesh(const FaceRange& faces, + TriangleMesh& tmesh, + const NamedParameters& np = parameters::default_values()) { - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::edge_descriptor edge_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - typedef typename GetGeomTraits::type GeomTraits; - typedef typename GetVertexPointMap::type VertexPointMap; - - // We need a default pmap that is not just 'constant_pmap(false)' because if an edge is constrained, - // its vertices are constrained. - typedef CGAL::dynamic_vertex_property_t Vertex_property_tag; - typedef typename boost::property_map::type Default_VCMap; - typedef typename internal_np::Lookup_named_param_def ::type VCMap; - - typedef typename internal_np::Lookup_named_param_def // default - > ::type ECMap; - - typedef internal::Area_smoother Area_optimizer; - typedef internal::Mesh_smoother Area_smoother; - typedef internal::Delaunay_edge_flipper Delaunay_flipper; - - typedef internal::Angle_smoother Angle_optimizer; - typedef internal::Mesh_smoother Angle_smoother; - - typedef typename GeomTraits::Triangle_3 Triangle; - typedef std::vector Triangle_container; - - typedef CGAL::AABB_triangle_primitive AABB_Primitive; - typedef CGAL::AABB_traits AABB_Traits; - typedef CGAL::AABB_tree Tree; - - if(std::begin(faces) == std::end(faces)) - return; - - using parameters::choose_parameter; - using parameters::get_parameter; - - // named parameters - GeomTraits gt = choose_parameter(get_parameter(np, internal_np::geom_traits)); - VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point), - get_property_map(CGAL::vertex_point, tmesh)); - - const bool use_angle_smoothing = choose_parameter(get_parameter(np, internal_np::use_angle_smoothing), true); - bool use_area_smoothing = choose_parameter(get_parameter(np, internal_np::use_area_smoothing), true); - -#ifndef CGAL_PMP_USE_CERES_SOLVER -#ifdef CGAL_PMP_SMOOTHING_DEBUG - if(use_area_smoothing) - std::cerr << "Requested area-based smoothing requires the Ceres Library, which is not available." << std::endl; -#endif - use_area_smoothing = false; -#endif - -#ifdef CGAL_PMP_SMOOTHING_DEBUG - if(!use_angle_smoothing && !use_area_smoothing) - std::cerr << "Called PMP::smooth_mesh() without any smoothing method selected or available" << std::endl; -#endif - - unsigned int nb_iterations = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 1); - const bool do_project = choose_parameter(get_parameter(np, internal_np::do_project), true); - const bool use_safety_constraints = choose_parameter(get_parameter(np, internal_np::use_safety_constraints), true); - const bool use_Delaunay_flips = choose_parameter(get_parameter(np, internal_np::use_Delaunay_flips), true); - - VCMap vcmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained), - get(Vertex_property_tag(), tmesh)); - - // If it's the default vcmap, manually set everything to false because the dynamic pmap has no default initialization - if((std::is_same::value)) - { - for(vertex_descriptor v : vertices(tmesh)) - put(vcmap, v, false); - } - - ECMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained), - Static_boolean_property_map()); - - // a constrained edge has constrained extremities - for(face_descriptor f : faces) - { - if(f == boost::graph_traits::null_face()) - continue; - - for(halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(f, tmesh), tmesh)) - { - if(get(ecmap, edge(h, tmesh))) - { - put(vcmap, source(h, tmesh), true); - put(vcmap, target(h, tmesh), true); - } - } - } - - // Construct the AABB tree (if needed for reprojection) - std::vector input_triangles; - - if(do_project) - { - input_triangles.reserve(faces.size()); - - for(face_descriptor f : faces) - { - halfedge_descriptor h = halfedge(f, tmesh); - if(is_border(h, tmesh)) // should not happen, but just in case - continue; - - input_triangles.push_back(gt.construct_triangle_3_object()(get(vpmap, source(h, tmesh)), - get(vpmap, target(h, tmesh)), - get(vpmap, target(next(h, tmesh), tmesh)))); - } - } - - Tree aabb_tree(input_triangles.begin(), input_triangles.end()); - - // Setup the working ranges and check some preconditions - Angle_smoother angle_smoother(tmesh, vpmap, vcmap, gt); - Area_smoother area_smoother(tmesh, vpmap, vcmap, gt); - Delaunay_flipper delaunay_flipper(tmesh, vpmap, ecmap, gt); - - if(use_angle_smoothing) - angle_smoother.init_smoothing(faces); - - if(use_area_smoothing) - area_smoother.init_smoothing(faces); - - for(unsigned int i=0; i -void smooth_mesh(TriangleMesh& tmesh, const CGAL_NP_CLASS& np = parameters::default_values()) +CGAL_DEPRECATED void smooth_mesh(TriangleMesh& tmesh, const CGAL_NP_CLASS& np = parameters::default_values()) { smooth_mesh(faces(tmesh), tmesh, np); } -template -void angles_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output) -{ - internal::Quality_evaluator evaluator(tmesh, traits); - evaluator.gather_angles(); - evaluator.extract_angles(output); -} - -template -void areas_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output) -{ - internal::Quality_evaluator evaluator(tmesh, traits); - evaluator.measure_areas(); - evaluator.extract_areas(output); -} - -template -void aspect_ratio_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output) -{ - internal::Quality_evaluator evaluator(tmesh, traits); - evaluator.calc_aspect_ratios(); - evaluator.extract_aspect_ratios(output); -} -///\endcond SKIP_IN_MANUAL } // namespace Polygon_mesh_processing } // namespace CGAL +#endif //#ifndef CGAL_NO_DEPRECATED_CODE + #endif // CGAL_POLYGON_MESH_PROCESSING_SMOOTH_MESH_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h new file mode 100644 index 00000000000..f4462e75057 --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/tangential_relaxation.h @@ -0,0 +1,326 @@ +// Copyright (c) 2015-2021 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) : Jane Tournois + + +#ifndef CGAL_POLYGON_MESH_PROCESSING_TANGENTIAL_RELAXATION_H +#define CGAL_POLYGON_MESH_PROCESSING_TANGENTIAL_RELAXATION_H + +#include + +#include +#include + +#include +#include + +#include +#include +#include + +namespace CGAL { +namespace Polygon_mesh_processing { + +namespace internal { +struct Allow_all_moves{ + template + constexpr inline bool operator()(vertex_descriptor, Point_3, Point_3) const + { + return true; + } +}; +} // internal namespace + + +/*! +* \ingroup PMP_meshing_grp +* applies an iterative area-based tangential smoothing to the given range of vertices. +* Each vertex `v` is relocated to its gravity-weighted centroid, and the relocation vector +* is projected back to the tangent plane to the surface at `v`, iteratively. +* The connectivity remains unchanged. +* +* @tparam TriangleMesh model of `FaceGraph` and `VertexListGraph`. +* The descriptor types `boost::graph_traits::%face_descriptor` +* and `boost::graph_traits::%halfedge_descriptor` must be +* models of `Hashable`. +* @tparam VertexRange range of `boost::graph_traits::%vertex_descriptor`, +* model of `Range`. Its iterator type is `ForwardIterator`. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" +* +* @param vertices the range of vertices which will be relocated by relaxation +* @param tm the triangle mesh to which `vertices` belong +* @param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below +* +* \cgalNamedParamsBegin +* \cgalParamNBegin{vertex_point_map} +* \cgalParamDescription{a property map associating points to the vertices of `tm`} +* \cgalParamType{a class model of `ReadWritePropertyMap` 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{geom_traits} +* \cgalParamDescription{an instance of a geometric traits class} +* \cgalParamType{a class model of `Kernel`} +* \cgalParamDefault{a \cgal Kernel deduced from the `Point_3` type, using `CGAL::Kernel_traits`} +* \cgalParamExtra{The geometric traits class must be compatible with the vertex `Point_3` type.} +* \cgalParamExtra{Exact constructions kernels are not supported by this function.} +* \cgalParamNEnd +* +* \cgalParamNBegin{number_of_iterations} +* \cgalParamDescription{the number of smoothing iterations} +* \cgalParamType{unsigned int} +* \cgalParamDefault{`1`} +* \cgalParamNEnd +* +* \cgalParamNBegin{edge_is_constrained_map} +* \cgalParamDescription{a property map containing the constrained-or-not status of each edge of `tm`. +* The endpoints of a constrained edge cannot be moved by relaxation.} +* \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits::%edge_descriptor` +* as key type and `bool` as value type. It must be default constructible.} +* \cgalParamDefault{a default property map where no edges are constrained} +* \cgalParamExtra{Boundary edges are always considered as constrained edges.} +* \cgalParamNEnd +* +* \cgalParamNBegin{vertex_is_constrained_map} +* \cgalParamDescription{a property map containing the constrained-or-not status of each vertex of `tm`. +* A constrained vertex cannot be modified during relaxation.} +* \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` +* as key type and `bool` as value type. It must be default constructible.} +* \cgalParamDefault{a default property map where no vertices are constrained} +* \cgalParamNEnd +* +* \cgalParamNBegin{relax_constraints} +* \cgalParamDescription{If `true`, the end vertices of the edges set as constrained +* in `edge_is_constrained_map` and boundary edges move along the +* constrained polylines they belong to.} +* \cgalParamType{Boolean} +* \cgalParamDefault{`false`} +* \cgalParamNEnd +* +* \cgalParamNBegin{allow_move_functor} +* \cgalParamDescription{A function object used to determinate if a vertex move should be allowed or not} +* \cgalParamType{Unary functor that provides `bool operator()(vertex_descriptor v, Point_3 src, Point_3 tgt)` returning `true` +* if the vertex `v` can be moved from `src` to `tgt`; %Point_3` being the value type of the vertex point map } +* \cgalParamDefault{If not provided, all moves are allowed.} +* \cgalParamNEnd +* +* \cgalNamedParamsEnd +* +* \todo check if it should really be a triangle mesh or if a polygon mesh is fine +*/ +template +void tangential_relaxation(const VertexRange& vertices, + TriangleMesh& tm, + const NamedParameters& np = parameters::default_values()) +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + + using parameters::get_parameter; + using parameters::choose_parameter; + + typedef typename GetGeomTraits::type GT; + + typedef typename GetVertexPointMap::type VPMap; + VPMap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), + get_property_map(vertex_point, tm)); + + typedef Static_boolean_property_map Default_ECM; + typedef typename internal_np::Lookup_named_param_def < + internal_np::edge_is_constrained_t, + NamedParameters, + Static_boolean_property_map // default (no constraint) + > ::type ECM; + ECM ecm = choose_parameter(get_parameter(np, internal_np::edge_is_constrained), + Default_ECM()); + + typedef typename internal_np::Lookup_named_param_def < + internal_np::vertex_is_constrained_t, + NamedParameters, + Static_boolean_property_map // default (no constraint) + > ::type VCM; + VCM vcm = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained), + Static_boolean_property_map()); + + const bool relax_constraints = choose_parameter(get_parameter(np, internal_np::relax_constraints), false); + const unsigned int nb_iterations = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 1); + + typedef typename GT::Vector_3 Vector_3; + typedef typename GT::Point_3 Point_3; + + auto check_normals = [&](vertex_descriptor v) + { + bool first_run = true; + Vector_3 prev = NULL_VECTOR, first = NULL_VECTOR; + halfedge_descriptor first_h = boost::graph_traits::null_halfedge(); + for (halfedge_descriptor hd : CGAL::halfedges_around_target(v, tm)) + { + if (is_border(hd, tm)) continue; + + Vector_3 n = compute_face_normal(face(hd, tm), tm, np); + if (n == CGAL::NULL_VECTOR) //for degenerate faces + continue; + + if (first_run) + { + first_run = false; + first = n; + first_h = hd; + } + else + { + if (!get(ecm, edge(hd, tm))) + if (to_double(n * prev) <= 0) + return false; + } + prev = n; + } + if (!get(ecm, edge(first_h, tm))) + if (to_double(first * prev) <= 0) + return false; + + return true; + }; + + typedef typename internal_np::Lookup_named_param_def < + internal_np::allow_move_functor_t, + NamedParameters, + internal::Allow_all_moves// default + > ::type Shall_move; + Shall_move shall_move = choose_parameter(get_parameter(np, internal_np::allow_move_functor), + internal::Allow_all_moves()); + + for (unsigned int nit = 0; nit < nb_iterations; ++nit) + { +#ifdef CGAL_PMP_TANGENTIAL_RELAXATION_VERBOSE + std::cout << "\r\t(Tangential relaxation iteration " << (nit + 1) << " / "; + std::cout << nb_iterations << ") "; + std::cout.flush(); +#endif + + typedef std::tuple VNP; + std::vector< VNP > barycenters; + + // at each vertex, compute vertex normal + std::unordered_map vnormals; + compute_vertex_normals(tm, boost::make_assoc_property_map(vnormals), np); + + // at each vertex, compute barycenter of neighbors + for(vertex_descriptor v : vertices) + { + if (get(vcm, v) || CGAL::internal::is_isolated(v, tm)) + continue; + + // collect hedges to detect if we have to handle boundary cases + std::vector interior_hedges, border_halfedges; + for(halfedge_descriptor h : halfedges_around_target(v, tm)) + { + if (is_border_edge(h, tm) || get(ecm, edge(h, tm))) + border_halfedges.push_back(h); + else + interior_hedges.push_back(h); + } + + if (border_halfedges.empty()) + { + const Vector_3& vn = vnormals.at(v); + Vector_3 move = CGAL::NULL_VECTOR; + unsigned int star_size = 0; + for(halfedge_descriptor h :interior_hedges) + { + move = move + Vector_3(get(vpm, v), get(vpm, source(h, tm))); + ++star_size; + } + CGAL_assertion(star_size > 0); //isolated vertices have already been discarded + move = (1. / static_cast(star_size)) * move; + + barycenters.emplace_back(v, vn, get(vpm, v) + move); + } + else + { + if (!relax_constraints) continue; + Vector_3 vn(NULL_VECTOR); + + if (border_halfedges.size() == 2)// corners are constrained + { + vertex_descriptor ph0 = source(border_halfedges[0], tm); + vertex_descriptor ph1 = source(border_halfedges[1], tm); + double dot = to_double(Vector_3(get(vpm, v), get(vpm, ph0)) + * Vector_3(get(vpm, v), get(vpm, ph1))); + // \todo shouldn't it be an input parameter? + //check squared cosine is < 0.25 (~120 degrees) + if (0.25 < dot*dot / ( squared_distance(get(vpm,ph0), get(vpm, v)) * + squared_distance(get(vpm,ph1), get(vpm, v))) ) + barycenters.emplace_back(v, vn, barycenter(get(vpm, ph0), 0.25, get(vpm, ph1), 0.25, get(vpm, v), 0.5)); + } + } + } + + // compute moves + typedef std::pair VP_pair; + std::vector< std::pair > new_locations; + new_locations.reserve(barycenters.size()); + for(const VNP& vnp : barycenters) + { + vertex_descriptor v = std::get<0>(vnp); + const Point_3& pv = get(vpm, v); + const Vector_3& nv = std::get<1>(vnp); + const Point_3& qv = std::get<2>(vnp); //barycenter at v + + new_locations.emplace_back(v, qv + (nv * Vector_3(qv, pv)) * nv); + } + + // perform moves + for(const VP_pair& vp : new_locations) + { + const Point_3 initial_pos = get(vpm, vp.first); // make a copy on purpose + const Vector_3 move(initial_pos, vp.second); + + put(vpm, vp.first, vp.second); + + //check that no inversion happened + double frac = 1.; + while (frac > 0.03 //5 attempts maximum + && ( !check_normals(vp.first) + || !shall_move(vp.first, initial_pos, get(vpm, vp.first)))) //if a face has been inverted + { + frac = 0.5 * frac; + put(vpm, vp.first, initial_pos + frac * move);//shorten the move by 2 + } + if (frac <= 0.02) + put(vpm, vp.first, initial_pos);//cancel move + } + }//end for loop (nit == nb_iterations) + +#ifdef CGAL_PMP_TANGENTIAL_RELAXATION_VERBOSE + std::cout << "\rTangential relaxation : " + << nb_iterations << " iterations done." << std::endl; +#endif +} + +/*! +* \ingroup PMP_meshing_grp +* applies `tangential_relaxation()` to all the vertices of `tm`. +*/ +template +void tangential_relaxation(TriangleMesh& tm, const CGAL_NP_CLASS& np = parameters::default_values()) +{ + tangential_relaxation(vertices(tm), tm, np); +} + +} } // CGAL::Polygon_mesh_processing + +#endif //CGAL_POLYGON_MESH_PROCESSING_TANGENTIAL_RELAXATION_H diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_smoothing.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_smoothing.cpp index 238c17768b4..986bdca5b53 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_smoothing.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_smoothing.cpp @@ -2,7 +2,8 @@ #include #include -#include +#include +#include #include @@ -36,8 +37,8 @@ void test_smoothing(const std::string filename) Mesh mesh; read_mesh(filename, mesh); - PMP::smooth_mesh(mesh); - PMP::smooth_mesh(mesh, CGAL::parameters::number_of_iterations(10)); + PMP::angle_and_area_smoothing(mesh); + PMP::angle_and_area_smoothing(mesh, CGAL::parameters::number_of_iterations(10)); } template @@ -46,9 +47,9 @@ void test_angle_smoothing(const std::string filename) Mesh mesh; read_mesh(filename, mesh); - PMP::smooth_mesh(mesh); - PMP::smooth_mesh(mesh, CGAL::parameters::number_of_iterations(10) - .use_area_smoothing(false)); + PMP::angle_and_area_smoothing(mesh); + PMP::angle_and_area_smoothing(mesh, CGAL::parameters::number_of_iterations(10) + .use_area_smoothing(false)); } template @@ -57,9 +58,9 @@ void test_area_smoothing(const std::string filename) Mesh mesh; read_mesh(filename, mesh); - PMP::smooth_mesh(mesh); - PMP::smooth_mesh(mesh, CGAL::parameters::number_of_iterations(10) - .use_angle_smoothing(false)); + PMP::angle_and_area_smoothing(mesh); + PMP::angle_and_area_smoothing(mesh, CGAL::parameters::number_of_iterations(10) + .use_angle_smoothing(false)); } template @@ -68,8 +69,8 @@ void test_angle_smoothing_without_projection(const std::string filename) Mesh mesh; read_mesh(filename, mesh); - PMP::smooth_mesh(mesh, CGAL::parameters::do_project(false) - .use_area_smoothing(false)); + PMP::angle_and_area_smoothing(mesh, CGAL::parameters::do_project(false) + .use_area_smoothing(false)); } template @@ -78,8 +79,38 @@ void test_area_smoothing_without_projection(const std::string filename) Mesh mesh; read_mesh(filename, mesh); - PMP::smooth_mesh(mesh, CGAL::parameters::do_project(false) - .use_angle_smoothing(false)); + PMP::angle_and_area_smoothing(mesh, CGAL::parameters::do_project(false) + .use_angle_smoothing(false)); +} + +template +void test_tangential_relaxation(const std::string filename) +{ + Mesh mesh; + read_mesh(filename, mesh); + PMP::tangential_relaxation( + vertices(mesh), + mesh, + CGAL::parameters::number_of_iterations(4) + .relax_constraints(false)); + PMP::tangential_relaxation(vertices(mesh), mesh); + PMP::tangential_relaxation(mesh); + + // test allow_move_functor + typename boost::property_map >::type vpm = + get(CGAL::dynamic_vertex_property_t(), mesh); + for (auto v : vertices(mesh)) + put(vpm, v, get(CGAL::vertex_point, mesh, v)); + auto no_move = [](typename boost::graph_traits::vertex_descriptor, Point, Point) + { + return false; + }; + PMP::tangential_relaxation(vertices(mesh), mesh, CGAL::parameters::allow_move_functor(no_move) + .vertex_point_map(vpm)); + for (auto v : vertices(mesh)) + { + assert(get(vpm, v) == get(CGAL::vertex_point, mesh, v)); + } } template @@ -104,7 +135,7 @@ void test_constrained_vertices(const std::string filename) CGAL::Boolean_property_map > vcmap(selected_vertices); - PMP::smooth_mesh(mesh, CGAL::parameters::vertex_is_constrained_map(vcmap)); + PMP::angle_and_area_smoothing(mesh, CGAL::parameters::vertex_is_constrained_map(vcmap)); for(vertex_descriptor v : vertices(mesh)) { @@ -127,6 +158,7 @@ int main(int /*argc*/, char** /*argv*/) test_angle_smoothing_without_projection(filename_elephant); test_area_smoothing_without_projection(filename_mannequin); test_constrained_vertices(filename_elephant); + test_tangential_relaxation(filename_elephant); // test with Polyhedron test_smoothing(filename_elephant); @@ -135,6 +167,7 @@ int main(int /*argc*/, char** /*argv*/) test_angle_smoothing_without_projection(filename_elephant); test_area_smoothing_without_projection(filename_mannequin); test_constrained_vertices(filename_mannequin); + test_tangential_relaxation(filename_elephant); std::cout << "Done!" << std::endl; return EXIT_SUCCESS; diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt b/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt index 26c24ce0875..76b60fbd535 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt @@ -38,7 +38,7 @@ if(TARGET CGAL::Eigen3_support) demo_framework CGAL::Eigen3_support) # The smoothing plugin can still do some things, even if Ceres is not found - qt5_wrap_ui( smoothingUI_FILES Smoothing_plugin.ui) + qt5_wrap_ui( smoothingUI_FILES Smoothing_plugin.ui Smoothing_tangential_relaxation.ui) polyhedron_demo_plugin(smoothing_plugin Smoothing_plugin ${smoothingUI_FILES}) target_link_libraries(smoothing_plugin PUBLIC scene_surface_mesh_item scene_selection_item CGAL::Eigen3_support) find_package(Ceres QUIET) diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Smoothing_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Smoothing_plugin.cpp index 15e41fe939b..8b08e92b62b 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Smoothing_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Smoothing_plugin.cpp @@ -8,14 +8,16 @@ #include #include -#include +#include #include +#include #include "Scene.h" #include "Scene_surface_mesh_item.h" #include "Scene_polyhedron_selection_item.h" #include "ui_Smoothing_plugin.h" +#include "ui_Smoothing_tangential_relaxation.h" #include #include @@ -69,11 +71,15 @@ public: connect(ui_widget.mesh_smoothing_button, SIGNAL(clicked()), this, SLOT(on_mesh_smoothing_clicked())); connect(ui_widget.shape_smoothing_button, SIGNAL(clicked()), this, SLOT(on_shape_smoothing_clicked())); + + actionRelax_ = new QAction(tr("Tangential Relaxation"), mw); + actionRelax_->setProperty("subMenuName", "Polygon Mesh Processing"); + connect(actionRelax_, SIGNAL(triggered()), this, SLOT(tangential_relaxation_action())); } QList actions() const { - return QList() << actionSmoothing_; + return QList() << actionSmoothing_ << actionRelax_; } bool applicable(QAction*) const @@ -127,6 +133,23 @@ public: } } + Ui::Tangential_relaxation_dialog + relaxation_dialog(QDialog* dialog) + { + Ui::Tangential_relaxation_dialog ui; + ui.setupUi(dialog); + connect(ui.buttonBox, SIGNAL(accepted()), dialog, SLOT(accept())); + connect(ui.buttonBox, SIGNAL(rejected()), dialog, SLOT(reject())); + + ui.nbIterations_spinbox->setSingleStep(1); + ui.nbIterations_spinbox->setRange(1/*min*/, 1000/*max*/); + ui.nbIterations_spinbox->setValue(1); + + ui.smooth1D_checkbox->setChecked(true); + + return ui; + } + public Q_SLOTS: void smoothing_action() { @@ -145,6 +168,85 @@ public Q_SLOTS: } } + void tangential_relaxation_action() + { + const Scene_interface::Item_id index = scene->mainSelectionIndex(); + + Scene_facegraph_item* poly_item = + qobject_cast(scene->item(index)); + Scene_polyhedron_selection_item* selection_item = + qobject_cast(scene->item(index)); + + if (poly_item || selection_item) + { + if (selection_item && selection_item->selected_facets.empty()) + { + QMessageBox::warning(mw, "Empty Facets", "There are no selected facets. Aborting."); + return; + } + // Create dialog box + QDialog dialog(mw); + Ui::Tangential_relaxation_dialog ui = relaxation_dialog(&dialog); + + // Get values + int i = dialog.exec(); + if (i == QDialog::Rejected) + { + std::cout << "Tangential relaxation aborted" << std::endl; + return; + } + + unsigned int nb_iter = ui.nbIterations_spinbox->value(); + bool smooth_features = ui.smooth1D_checkbox->isChecked(); + + // wait cursor + QApplication::setOverrideCursor(Qt::WaitCursor); + QElapsedTimer time; + time.start(); + + FaceGraph& pmesh = (poly_item != nullptr) + ? *poly_item->polyhedron() + : *selection_item->polyhedron(); + + if (selection_item) + { + boost::unordered_set vset; + for (face_descriptor f : selection_item->selected_facets) + { + for(vertex_descriptor fv : CGAL::vertices_around_face(halfedge(f, pmesh), pmesh)) + vset.insert(fv); + } + + CGAL::Polygon_mesh_processing::tangential_relaxation( + vset, + pmesh, + CGAL::Polygon_mesh_processing::parameters::number_of_iterations(nb_iter) + .edge_is_constrained_map(selection_item->constrained_edges_pmap()) + .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()) + .relax_constraints(smooth_features)); + selection_item->polyhedron_item()->invalidateOpenGLBuffers(); + Q_EMIT selection_item->polyhedron_item()->itemChanged(); + selection_item->invalidateOpenGLBuffers(); + selection_item->setKeepSelectionValid(Scene_polyhedron_selection_item::None); + } + else if (poly_item) + { + CGAL::Polygon_mesh_processing::tangential_relaxation( + vertices(pmesh), + pmesh, + CGAL::Polygon_mesh_processing::parameters::number_of_iterations(nb_iter)); + + poly_item->invalidateOpenGLBuffers(); + Q_EMIT poly_item->itemChanged(); + } + + std::cout << "ok (" << time.elapsed() << " ms)" << std::endl; + } + + // default cursor + QApplication::restoreOverrideCursor(); + } + void on_mesh_smoothing_clicked() { const Scene_interface::Item_id index = scene->mainSelectionIndex(); @@ -175,7 +277,7 @@ public Q_SLOTS: if(poly_item) { - smooth_mesh(pmesh, parameters::do_project(projection) + angle_and_area_smoothing(pmesh, parameters::do_project(projection) .number_of_iterations(nb_iter) .vertex_is_constrained_map(vcmap) .use_safety_constraints(use_safety_measures) @@ -193,7 +295,7 @@ public Q_SLOTS: // No faces selected --> use all faces if(std::begin(selection_item->selected_facets) == std::end(selection_item->selected_facets)) { - smooth_mesh(pmesh, parameters::do_project(projection) + angle_and_area_smoothing(pmesh, parameters::do_project(projection) .number_of_iterations(nb_iter) .vertex_is_constrained_map(vcmap) .edge_is_constrained_map(selection_item->constrained_edges_pmap()) @@ -204,7 +306,7 @@ public Q_SLOTS: } else // some faces exist in the selection { - smooth_mesh(selection_item->selected_facets, pmesh, parameters::do_project(projection) + angle_and_area_smoothing(selection_item->selected_facets, pmesh, parameters::do_project(projection) .number_of_iterations(nb_iter) .vertex_is_constrained_map(vcmap) .edge_is_constrained_map(selection_item->constrained_edges_pmap()) @@ -284,6 +386,7 @@ public Q_SLOTS: private: QAction* actionSmoothing_; + QAction* actionRelax_; QDockWidget* dock_widget; Ui::Smoothing ui_widget; }; diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Smoothing_tangential_relaxation.ui b/Polyhedron/demo/Polyhedron/Plugins/PMP/Smoothing_tangential_relaxation.ui new file mode 100644 index 00000000000..fd9c4aefe1d --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Smoothing_tangential_relaxation.ui @@ -0,0 +1,177 @@ + + + Tangential_relaxation_dialog + + + true + + + + 0 + 0 + 376 + 206 + + + + Isotropic remeshing criteria + + + + + + + 15 + 75 + true + + + + NO OBJECT + + + + + + + No size + + + + + + + Tangential Relaxation + + + + + + Number of iterations + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + nbIterations_spinbox + + + + + + + + 110 + 0 + + + + 1 + + + + + + + Allow 1D relaxation + (along borders/constraints) + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + + Qt::Vertical + + + QSizePolicy::Maximum + + + + 20 + 40 + + + + + + + + + + + + + + + + + + Qt::Vertical + + + QSizePolicy::MinimumExpanding + + + + 0 + 0 + + + + + + + + Qt::Horizontal + + + QDialogButtonBox::Cancel|QDialogButtonBox::Ok + + + + + + + nbIterations_spinbox + smooth1D_checkbox + + + + + buttonBox + accepted() + Tangential_relaxation_dialog + accept() + + + 397 + 333 + + + 157 + 195 + + + + + buttonBox + rejected() + Tangential_relaxation_dialog + reject() + + + 397 + 333 + + + 286 + 195 + + + + + 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 d597cf0b5d5..d9fe8962785 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -90,6 +90,7 @@ CGAL_add_named_parameter(overlap_test_t, overlap_test, do_overlap_test_of_bounde CGAL_add_named_parameter(preserve_genus_t, preserve_genus, preserve_genus) CGAL_add_named_parameter(apply_per_connected_component_t, apply_per_connected_component, apply_per_connected_component) CGAL_add_named_parameter(projection_functor_t, projection_functor, projection_functor) +CGAL_add_named_parameter(allow_move_functor_t, allow_move_functor, allow_move_functor) CGAL_add_named_parameter(throw_on_self_intersection_t, throw_on_self_intersection, throw_on_self_intersection) CGAL_add_named_parameter(clip_volume_t, clip_volume, clip_volume) CGAL_add_named_parameter(use_compact_clipper_t, use_compact_clipper, use_compact_clipper)