Merge remote-tracking branch 'sloriot/PMP-relax'

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<Point_3>& poly) -> std::size_t
                        { return size + poly.size(); });
 

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)
 -----------
 

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()`
 

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 <i>blobby</i>, 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.

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

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")

Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_smoothing_example.cpp

@@ -1,7 +1,7 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Surface_mesh.h>
 
-#include <CGAL/Polygon_mesh_processing/smooth_mesh.h>
+#include <CGAL/Polygon_mesh_processing/angle_and_area_smoothing.h>
 #include <CGAL/Polygon_mesh_processing/detect_features.h>
 #include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
 
@@ -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)
+  PMP::angle_and_area_smoothing(mesh, CGAL::parameters::number_of_iterations(nb_iterations)
-                                          .use_safety_constraints(false) // authorize all moves
+                                                       .use_safety_constraints(false) // authorize all moves
-                                          .edge_is_constrained_map(eif));
+                                                       .edge_is_constrained_map(eif));
 
   CGAL::IO::write_polygon_mesh("mesh_smoothed.off", mesh, CGAL::parameters::stream_precision(17));
 

Polygon_mesh_processing/examples/Polygon_mesh_processing/tangential_relaxation_example.cpp

@@ -0,0 +1,35 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/tangential_relaxation.h>
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+
+#include <fstream>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel   K;
+typedef CGAL::Surface_mesh<K::Point_3>                        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;
+}

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 <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
+
+#include <CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h>
+#include <CGAL/Polygon_mesh_processing/internal/Smoothing/smoothing_evaluation.h>
+#include <CGAL/Polygon_mesh_processing/self_intersections.h>
+
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#include <CGAL/property_map.h>
+
+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<TriangleMesh>::%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<TriangleMesh>::%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<TriangleMesh>::%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<TriangleMesh>::%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<typename TriangleMesh, typename FaceRange, typename NamedParameters = parameters::Default_named_parameters>
+void angle_and_area_smoothing(const FaceRange& faces,
+                              TriangleMesh& tmesh,
+                              const NamedParameters& np = parameters::default_values())
+{
+  typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor               vertex_descriptor;
+  typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor             halfedge_descriptor;
+  typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor                 edge_descriptor;
+  typedef typename boost::graph_traits<TriangleMesh>::face_descriptor                 face_descriptor;
+
+  typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type                 GeomTraits;
+  typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::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<bool>                                       Vertex_property_tag;
+  typedef typename boost::property_map<TriangleMesh, Vertex_property_tag>::type       Default_VCMap;
+  typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_is_constrained_t,
+                                                 NamedParameters,
+                                                 Default_VCMap
+                                                 > ::type                             VCMap;
+
+  typedef typename internal_np::Lookup_named_param_def<internal_np::edge_is_constrained_t,
+                                                 NamedParameters,
+                                                 Static_boolean_property_map<edge_descriptor, false> // default
+                                                 > ::type                             ECMap;
+
+  typedef internal::Area_smoother<TriangleMesh, VertexPointMap, GeomTraits>           Area_optimizer;
+  typedef internal::Mesh_smoother<Area_optimizer, TriangleMesh,
+                                  VertexPointMap, VCMap, GeomTraits>                  Area_smoother;
+  typedef internal::Delaunay_edge_flipper<TriangleMesh, VertexPointMap,
+                                          ECMap, GeomTraits>                          Delaunay_flipper;
+
+  typedef internal::Angle_smoother<TriangleMesh, VertexPointMap, GeomTraits>          Angle_optimizer;
+  typedef internal::Mesh_smoother<Angle_optimizer, TriangleMesh,
+                                  VertexPointMap, VCMap, GeomTraits>                  Angle_smoother;
+
+  typedef typename GeomTraits::Triangle_3                                             Triangle;
+  typedef std::vector<Triangle>                                                       Triangle_container;
+
+  typedef CGAL::AABB_triangle_primitive<GeomTraits,
+                                        typename Triangle_container::iterator>        AABB_Primitive;
+  typedef CGAL::AABB_traits<GeomTraits, AABB_Primitive>                               AABB_Traits;
+  typedef CGAL::AABB_tree<AABB_Traits>                                                Tree;
+
+  if(std::begin(faces) == std::end(faces))
+    return;
+
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  // named parameters
+  GeomTraits gt = choose_parameter<GeomTraits>(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<VCMap, Default_VCMap>::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<edge_descriptor, false>());
+
+  // a constrained edge has constrained extremities
+  for(face_descriptor f : faces)
+  {
+    if(f == boost::graph_traits<TriangleMesh>::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<Triangle> 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<nb_iterations; ++i)
+  {
+#ifdef CGAL_PMP_SMOOTHING_DEBUG
+    std::cout << "Iteration #" << i << std::endl;
+#endif
+
+    if(use_area_smoothing)
+    {
+#ifdef CGAL_PMP_SMOOTHING_DEBUG
+      std::cout << "Smooth areas..." << std::endl;
+#endif
+
+      // First apply area smoothing...
+      area_smoother.optimize(use_safety_constraints /*check for bad faces*/,
+                             false /*apply moves as soon as they're calculated*/,
+                             false /*do not enforce a minimum angle improvement*/);
+      if(do_project)
+      {
+        if(use_safety_constraints && does_self_intersect(tmesh))
+        {
+#ifdef CGAL_PMP_SMOOTHING_DEBUG
+          std::cerr << "Cannot re-project as there are self-intersections in the mesh!\n";
+#endif
+          break;
+        }
+
+        area_smoother.project_to_surface(aabb_tree);
+      }
+
+      if(use_Delaunay_flips)
+        delaunay_flipper(faces);
+    }
+
+    // ... then angle smoothing
+    if(use_angle_smoothing)
+    {
+#ifdef CGAL_PMP_SMOOTHING_DEBUG
+      std::cout << "Smooth angles..." << std::endl;
+#endif
+
+      angle_smoother.optimize(use_safety_constraints /*check for bad faces*/,
+                              true /*apply all moves at once*/,
+                              use_safety_constraints /*check if the min angle is improved*/);
+
+      if(do_project)
+      {
+        if(use_safety_constraints && does_self_intersect(tmesh))
+        {
+#ifdef CGAL_PMP_SMOOTHING_DEBUG
+          std::cerr << "Can't do re-projection, there are self-intersections in the mesh!\n";
+#endif
+          break;
+        }
+
+        angle_smoother.project_to_surface(aabb_tree);
+      }
+    }
+  }
+}
+
+///\cond SKIP_IN_MANUAL
+template <typename TriangleMesh, typename CGAL_NP_TEMPLATE_PARAMETERS>
+void angle_and_area_smoothing(TriangleMesh& tmesh, const CGAL_NP_CLASS& np = parameters::default_values())
+{
+  angle_and_area_smoothing(faces(tmesh), tmesh, np);
+}
+
+template<typename TriangleMesh, typename GeomTraits, typename Stream>
+void angles_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
+{
+  internal::Quality_evaluator<TriangleMesh, GeomTraits> evaluator(tmesh, traits);
+  evaluator.gather_angles();
+  evaluator.extract_angles(output);
+}
+
+template<typename TriangleMesh, typename GeomTraits, typename Stream>
+void areas_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
+{
+  internal::Quality_evaluator<TriangleMesh, GeomTraits> evaluator(tmesh, traits);
+  evaluator.measure_areas();
+  evaluator.extract_areas(output);
+}
+
+template<typename TriangleMesh, typename GeomTraits, typename Stream>
+void aspect_ratio_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
+{
+  internal::Quality_evaluator<TriangleMesh, GeomTraits> 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

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h

@@ -22,6 +22,7 @@
 #include <CGAL/Polygon_mesh_processing/measure.h>
 #include <CGAL/Polygon_mesh_processing/connected_components.h>
 #include <CGAL/Polygon_mesh_processing/shape_predicates.h>
+#include <CGAL/Polygon_mesh_processing/tangential_relaxation.h>
 
 #include <CGAL/AABB_tree.h>
 #include <CGAL/AABB_traits.h>
@@ -43,6 +44,7 @@
 #include <memory>
 #include <boost/container/flat_set.hpp>
 #include <boost/optional.hpp>
+#include <boost/property_map/function_property_map.hpp>
 
 #include <map>
 #include <list>
@@ -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*/
+    void tangential_relaxation_impl(const bool relax_constraints/*1d smoothing*/
-                             , const unsigned int nb_iterations)
+                                   , 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
+        return this->is_constrained(e);
-        std::cout << "\r\t(iteration " << (nit + 1) << " / ";
+      };
-        std::cout << nb_iterations << ") ";
+      auto constrained_edges_pmap
-        std::cout.flush();
+        = boost::make_function_property_map<edge_descriptor>(edge_constraint);
-#endif
+
-        typedef std::tuple<vertex_descriptor, Vector_3, Point> VNP;
+      // property map of constrained vertices for relaxation
-        std::vector< VNP > barycenters;
+      auto vertex_constraint = [&](const vertex_descriptor v)
-        // at each vertex, compute vertex normal
+      {
-        // at each vertex, compute barycenter of neighbors
+        for (halfedge_descriptor h : halfedges_around_target(v, mesh_))
-        for(vertex_descriptor v : vertices(mesh_))
         {
-          if (is_constrained(v) || is_isolated(v))
+          Halfedge_status s = status(h);
-            continue;
+          if ( s == PATCH
-
+            || s == PATCH_BORDER
-          else if (is_on_patch(v))
+            || status(opposite(h, mesh_)) == PATCH_BORDER)
-          {
+            return false;
-            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<halfedge_descriptor> 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]))) );
-            }
-          }
         }
+        return true;
+      };
+      auto constrained_vertices_pmap
+        = boost::make_function_property_map<vertex_descriptor>(vertex_constraint);
 
-        // compute moves
+      tangential_relaxation(
-        typedef std::pair<vertex_descriptor, Point> VP_pair;
+        vertices(mesh_),
-        std::vector< std::pair<vertex_descriptor, Point> > new_locations;
+        mesh_,
-        new_locations.reserve(barycenters.size());
+        CGAL::parameters::number_of_iterations(nb_iterations)
-        for(const VNP& vnp : barycenters)
+                         .vertex_point_map(vpmap_)
-        {
+                         .geom_traits(gt_)
-          vertex_descriptor v = std::get<0>(vnp);
+                         .edge_is_constrained_map(constrained_edges_pmap)
-          Point pv = get(vpmap_, v);
+                         .vertex_is_constrained_map(constrained_vertices_pmap)
-          const Vector_3& nv = std::get<1>(vnp);
+                         .relax_constraints(relax_constraints)
-          const Point& qv = std::get<2>(vnp); //barycenter at v
+      );
 
-          new_locations.push_back( std::make_pair(v, qv + (nv * Vector_3(qv, pv)) * nv) );
+      CGAL_assertion(!input_mesh_is_valid_ || is_valid_polygon_mesh(mesh_));
-        }
-
-        // 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)
 
 #ifdef CGAL_PMP_REMESHING_DEBUG
       debug_self_intersections();

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

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h

@@ -22,7 +22,7 @@
 #include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
 #include <CGAL/Polygon_mesh_processing/remesh.h>
 #include <CGAL/Polygon_mesh_processing/self_intersections.h>
-#include <CGAL/Polygon_mesh_processing/smooth_mesh.h>
+#include <CGAL/Polygon_mesh_processing/angle_and_area_smoothing.h>
 #include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
 #ifndef CGAL_PMP_REMOVE_SELF_INTERSECTION_NO_POLYHEDRAL_ENVELOPE_CHECK
 #include <CGAL/Polyhedral_envelope.h>
@@ -546,9 +546,11 @@ bool remove_self_intersections_with_smoothing(std::set<typename boost::graph_tra
                   constrain_sharp_edges, dihedral_angle, weak_DA, eif, local_vpm, gt);
 
   // @todo choice of number of iterations? Till convergence && max of 100?
-  Polygon_mesh_processing::smooth_mesh(faces(local_mesh), local_mesh, CP::edge_is_constrained_map(eif)
+  Polygon_mesh_processing::angle_and_area_smoothing(faces(local_mesh),
-                                                                         .number_of_iterations(100)
+                                                    local_mesh,
-                                                                         .use_safety_constraints(false));
+                                                    CP::edge_is_constrained_map(eif)
+                                                    .number_of_iterations(100)
+                                                    .use_safety_constraints(false));
 
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT
   CGAL::IO::write_polygon_mesh("results/post_smoothing_local_mesh.off", local_mesh, CGAL::parameters::stream_precision(17));

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h

@@ -16,352 +16,42 @@
 
 #include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
 
-#include <CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h>
+#define CGAL_DEPRECATED_HEADER "<CGAL/Polygon_mesh_processing/smooth_mesh.h>"
-#include <CGAL/Polygon_mesh_processing/internal/Smoothing/smoothing_evaluation.h>
+#define CGAL_REPLACEMENT_HEADER "<CGAL/Polygon_mesh_processing/angle_and_area_smoothing.h>"
-#include <CGAL/Polygon_mesh_processing/self_intersections.h>
+#include <CGAL/Installation/internal/deprecation_warning.h>
 
-#include <CGAL/Named_function_parameters.h>
+#include <CGAL/Polygon_mesh_processing/angle_and_area_smoothing.h>
-#include <CGAL/boost/graph/named_params_helper.h>
 
-#include <CGAL/property_map.h>
+#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.
+* \deprecated This function is deprecated since \cgal 5.5,
-*
+* `CGAL::angle_and_area_smoothing()` should be used instead.
-* 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<TriangleMesh>::%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<TriangleMesh>::%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<TriangleMesh>::%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<TriangleMesh>::%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()`
 */
 template<typename TriangleMesh, typename FaceRange, typename NamedParameters = parameters::Default_named_parameters>
-void smooth_mesh(const FaceRange& faces,
+CGAL_DEPRECATED void smooth_mesh(const FaceRange& faces,
-                 TriangleMesh& tmesh,
+                                 TriangleMesh& tmesh,
-                 const NamedParameters& np = parameters::default_values())
+                                 const NamedParameters& np = parameters::default_values())
 {
-  typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor               vertex_descriptor;
+  angle_and_area_smoothing(faces, tmesh, np);
-  typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor             halfedge_descriptor;
-  typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor                 edge_descriptor;
-  typedef typename boost::graph_traits<TriangleMesh>::face_descriptor                 face_descriptor;
-
-  typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type                 GeomTraits;
-  typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::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<bool>                                       Vertex_property_tag;
-  typedef typename boost::property_map<TriangleMesh, Vertex_property_tag>::type       Default_VCMap;
-  typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_is_constrained_t,
-                                                 NamedParameters,
-                                                 Default_VCMap
-                                                 > ::type                             VCMap;
-
-  typedef typename internal_np::Lookup_named_param_def<internal_np::edge_is_constrained_t,
-                                                 NamedParameters,
-                                                 Static_boolean_property_map<edge_descriptor, false> // default
-                                                 > ::type                             ECMap;
-
-  typedef internal::Area_smoother<TriangleMesh, VertexPointMap, GeomTraits>           Area_optimizer;
-  typedef internal::Mesh_smoother<Area_optimizer, TriangleMesh,
-                                  VertexPointMap, VCMap, GeomTraits>                  Area_smoother;
-  typedef internal::Delaunay_edge_flipper<TriangleMesh, VertexPointMap,
-                                          ECMap, GeomTraits>                          Delaunay_flipper;
-
-  typedef internal::Angle_smoother<TriangleMesh, VertexPointMap, GeomTraits>          Angle_optimizer;
-  typedef internal::Mesh_smoother<Angle_optimizer, TriangleMesh,
-                                  VertexPointMap, VCMap, GeomTraits>                  Angle_smoother;
-
-  typedef typename GeomTraits::Triangle_3                                             Triangle;
-  typedef std::vector<Triangle>                                                       Triangle_container;
-
-  typedef CGAL::AABB_triangle_primitive<GeomTraits,
-                                        typename Triangle_container::iterator>        AABB_Primitive;
-  typedef CGAL::AABB_traits<GeomTraits, AABB_Primitive>                               AABB_Traits;
-  typedef CGAL::AABB_tree<AABB_Traits>                                                Tree;
-
-  if(std::begin(faces) == std::end(faces))
-    return;
-
-  using parameters::choose_parameter;
-  using parameters::get_parameter;
-
-  // named parameters
-  GeomTraits gt = choose_parameter<GeomTraits>(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<VCMap, Default_VCMap>::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<edge_descriptor, false>());
-
-  // a constrained edge has constrained extremities
-  for(face_descriptor f : faces)
-  {
-    if(f == boost::graph_traits<TriangleMesh>::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<Triangle> 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<nb_iterations; ++i)
-  {
-#ifdef CGAL_PMP_SMOOTHING_DEBUG
-    std::cout << "Iteration #" << i << std::endl;
-#endif
-
-    if(use_area_smoothing)
-    {
-#ifdef CGAL_PMP_SMOOTHING_DEBUG
-      std::cout << "Smooth areas..." << std::endl;
-#endif
-
-      // First apply area smoothing...
-      area_smoother.optimize(use_safety_constraints /*check for bad faces*/,
-                             false /*apply moves as soon as they're calculated*/,
-                             false /*do not enforce a minimum angle improvement*/);
-      if(do_project)
-      {
-        if(use_safety_constraints && does_self_intersect(tmesh))
-        {
-#ifdef CGAL_PMP_SMOOTHING_DEBUG
-          std::cerr << "Cannot re-project as there are self-intersections in the mesh!\n";
-#endif
-          break;
-        }
-
-        area_smoother.project_to_surface(aabb_tree);
-      }
-
-      if(use_Delaunay_flips)
-        delaunay_flipper(faces);
-    }
-
-    // ... then angle smoothing
-    if(use_angle_smoothing)
-    {
-#ifdef CGAL_PMP_SMOOTHING_DEBUG
-      std::cout << "Smooth angles..." << std::endl;
-#endif
-
-      angle_smoother.optimize(use_safety_constraints /*check for bad faces*/,
-                              true /*apply all moves at once*/,
-                              use_safety_constraints /*check if the min angle is improved*/);
-
-      if(do_project)
-      {
-        if(use_safety_constraints && does_self_intersect(tmesh))
-        {
-#ifdef CGAL_PMP_SMOOTHING_DEBUG
-          std::cerr << "Can't do re-projection, there are self-intersections in the mesh!\n";
-#endif
-          break;
-        }
-
-        angle_smoother.project_to_surface(aabb_tree);
-      }
-    }
-  }
 }
 
 ///\cond SKIP_IN_MANUAL
 template <typename TriangleMesh, typename CGAL_NP_TEMPLATE_PARAMETERS>
-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<typename TriangleMesh, typename GeomTraits, typename Stream>
-void angles_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
-{
-  internal::Quality_evaluator<TriangleMesh, GeomTraits> evaluator(tmesh, traits);
-  evaluator.gather_angles();
-  evaluator.extract_angles(output);
-}
-
-template<typename TriangleMesh, typename GeomTraits, typename Stream>
-void areas_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
-{
-  internal::Quality_evaluator<TriangleMesh, GeomTraits> evaluator(tmesh, traits);
-  evaluator.measure_areas();
-  evaluator.extract_areas(output);
-}
-
-template<typename TriangleMesh, typename GeomTraits, typename Stream>
-void aspect_ratio_evaluation(TriangleMesh& tmesh, GeomTraits traits, Stream& output)
-{
-  internal::Quality_evaluator<TriangleMesh, GeomTraits> 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

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 <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
+
+#include <CGAL/Polygon_mesh_processing/compute_normal.h>
+#include <CGAL/property_map.h>
+
+#include <CGAL/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#include <tuple>
+#include <type_traits>
+#include <unordered_map>
+
+namespace CGAL {
+namespace Polygon_mesh_processing {
+
+namespace internal {
+struct Allow_all_moves{
+  template <class vertex_descriptor, class Point_3>
+  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<TriangleMesh>::%face_descriptor`
+*         and `boost::graph_traits<TriangleMesh>::%halfedge_descriptor` must be
+*         models of `Hashable`.
+* @tparam VertexRange range of `boost::graph_traits<TriangleMesh>::%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<TriangleMesh>::%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<TriangleMesh>::%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<PolygonMesh>::%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 <typename VertexRange, class TriangleMesh, class NamedParameters = parameters::Default_named_parameters>
+void tangential_relaxation(const VertexRange& vertices,
+                           TriangleMesh& tm,
+                           const NamedParameters& np = parameters::default_values())
+{
+  typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
+  typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
+  typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor edge_descriptor;
+
+  using parameters::get_parameter;
+  using parameters::choose_parameter;
+
+  typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type GT;
+
+  typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::type VPMap;
+  VPMap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
+                               get_property_map(vertex_point, tm));
+
+  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,
+      Static_boolean_property_map<edge_descriptor, false> // 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<vertex_descriptor, false> // default (no constraint)
+    > ::type VCM;
+  VCM vcm = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained),
+                             Static_boolean_property_map<vertex_descriptor, false>());
+
+  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<TriangleMesh>::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<vertex_descriptor, Vector_3, Point_3> VNP;
+    std::vector< VNP > barycenters;
+
+    // at each vertex, compute vertex normal
+    std::unordered_map<vertex_descriptor, Vector_3> 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<halfedge_descriptor> 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<double>(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<vertex_descriptor, Point_3> VP_pair;
+    std::vector< std::pair<vertex_descriptor, Point_3> > 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 <class TriangleMesh,
+          typename CGAL_NP_TEMPLATE_PARAMETERS>
+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

Polygon_mesh_processing/test/Polygon_mesh_processing/test_mesh_smoothing.cpp

@@ -2,7 +2,8 @@
 
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/Polyhedron_3.h>
-#include <CGAL/Polygon_mesh_processing/smooth_mesh.h>
+#include <CGAL/Polygon_mesh_processing/angle_and_area_smoothing.h>
+#include <CGAL/Polygon_mesh_processing/tangential_relaxation.h>
 
 #include <CGAL/property_map.h>
 
@@ -36,8 +37,8 @@ void test_smoothing(const std::string filename)
   Mesh mesh;
   read_mesh(filename, mesh);
 
-  PMP::smooth_mesh(mesh);
+  PMP::angle_and_area_smoothing(mesh);
-  PMP::smooth_mesh(mesh, CGAL::parameters::number_of_iterations(10));
+  PMP::angle_and_area_smoothing(mesh, CGAL::parameters::number_of_iterations(10));
 }
 
 template <typename Mesh>
@@ -46,9 +47,9 @@ void test_angle_smoothing(const std::string filename)
   Mesh mesh;
   read_mesh(filename, mesh);
 
-  PMP::smooth_mesh(mesh);
+  PMP::angle_and_area_smoothing(mesh);
-  PMP::smooth_mesh(mesh, CGAL::parameters::number_of_iterations(10)
+  PMP::angle_and_area_smoothing(mesh, CGAL::parameters::number_of_iterations(10)
-                                          .use_area_smoothing(false));
+                                                       .use_area_smoothing(false));
 }
 
 template <typename Mesh>
@@ -57,9 +58,9 @@ void test_area_smoothing(const std::string filename)
   Mesh mesh;
   read_mesh(filename, mesh);
 
-  PMP::smooth_mesh(mesh);
+  PMP::angle_and_area_smoothing(mesh);
-  PMP::smooth_mesh(mesh, CGAL::parameters::number_of_iterations(10)
+  PMP::angle_and_area_smoothing(mesh, CGAL::parameters::number_of_iterations(10)
-                                          .use_angle_smoothing(false));
+                                                       .use_angle_smoothing(false));
 }
 
 template <typename Mesh>
@@ -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)
+  PMP::angle_and_area_smoothing(mesh, CGAL::parameters::do_project(false)
-                                          .use_area_smoothing(false));
+                                                       .use_area_smoothing(false));
 }
 
 template <typename Mesh>
@@ -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)
+  PMP::angle_and_area_smoothing(mesh, CGAL::parameters::do_project(false)
-                                          .use_angle_smoothing(false));
+                                                       .use_angle_smoothing(false));
+}
+
+template<typename Mesh>
+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<Mesh, CGAL::dynamic_vertex_property_t<Point> >::type vpm =
+    get(CGAL::dynamic_vertex_property_t<Point>(), mesh);
+  for (auto v : vertices(mesh))
+    put(vpm, v, get(CGAL::vertex_point, mesh, v));
+  auto no_move = [](typename boost::graph_traits<Mesh>::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 <typename Mesh>
@@ -104,7 +135,7 @@ void test_constrained_vertices(const std::string filename)
 
   CGAL::Boolean_property_map<std::set<vertex_descriptor> > 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<SurfaceMesh>(filename_elephant);
   test_area_smoothing_without_projection<SurfaceMesh>(filename_mannequin);
   test_constrained_vertices<SurfaceMesh>(filename_elephant);
+  test_tangential_relaxation<SurfaceMesh>(filename_elephant);
 
   // test with Polyhedron
   test_smoothing<Polyhedron>(filename_elephant);
@@ -135,6 +167,7 @@ int main(int /*argc*/, char** /*argv*/)
   test_angle_smoothing_without_projection<Polyhedron>(filename_elephant);
   test_area_smoothing_without_projection<Polyhedron>(filename_mannequin);
   test_constrained_vertices<Polyhedron>(filename_mannequin);
+  test_tangential_relaxation<Polyhedron>(filename_elephant);
 
   std::cout << "Done!" << std::endl;
   return EXIT_SUCCESS;

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)

Polyhedron/demo/Polyhedron/Plugins/PMP/Smoothing_plugin.cpp

@@ -8,14 +8,16 @@
 #include <CGAL/property_map.h>
 #include <CGAL/Surface_mesh.h>
 
-#include <CGAL/Polygon_mesh_processing/smooth_mesh.h>
+#include <CGAL/Polygon_mesh_processing/angle_and_area_smoothing.h>
 #include <CGAL/Polygon_mesh_processing/smooth_shape.h>
+#include <CGAL/Polygon_mesh_processing/tangential_relaxation.h>
 
 #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 <QElapsedTimer>
 #include <QAction>
@@ -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<QAction*> actions() const
   {
-    return QList<QAction*>() << actionSmoothing_;
+    return QList<QAction*>() << 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_facegraph_item*>(scene->item(index));
+    Scene_polyhedron_selection_item* selection_item =
+      qobject_cast<Scene_polyhedron_selection_item*>(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<vertex_descriptor> 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;
 };

Polyhedron/demo/Polyhedron/Plugins/PMP/Smoothing_tangential_relaxation.ui

@@ -0,0 +1,177 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<ui version="4.0">
+ <class>Tangential_relaxation_dialog</class>
+ <widget class="QDialog" name="Tangential_relaxation_dialog">
+  <property name="enabled">
+   <bool>true</bool>
+  </property>
+  <property name="geometry">
+   <rect>
+    <x>0</x>
+    <y>0</y>
+    <width>376</width>
+    <height>206</height>
+   </rect>
+  </property>
+  <property name="windowTitle">
+   <string>Isotropic remeshing criteria</string>
+  </property>
+  <layout class="QVBoxLayout" name="verticalLayout">
+   <item>
+    <widget class="QLabel" name="objectName">
+     <property name="font">
+      <font>
+       <pointsize>15</pointsize>
+       <weight>75</weight>
+       <bold>true</bold>
+      </font>
+     </property>
+     <property name="text">
+      <string>NO OBJECT</string>
+     </property>
+    </widget>
+   </item>
+   <item>
+    <widget class="QLabel" name="objectNameSize">
+     <property name="text">
+      <string>No size</string>
+     </property>
+    </widget>
+   </item>
+   <item>
+    <widget class="QGroupBox" name="groupBox">
+     <property name="title">
+      <string>Tangential Relaxation</string>
+     </property>
+     <layout class="QGridLayout" name="gridLayout_2" rowstretch="0,0,0" columnstretch="0,0">
+      <item row="0" column="0">
+       <widget class="QLabel" name="nbIterations_label">
+        <property name="text">
+         <string>Number of iterations</string>
+        </property>
+        <property name="alignment">
+         <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
+        </property>
+        <property name="buddy">
+         <cstring>nbIterations_spinbox</cstring>
+        </property>
+       </widget>
+      </item>
+      <item row="0" column="1">
+       <widget class="QSpinBox" name="nbIterations_spinbox">
+        <property name="minimumSize">
+         <size>
+          <width>110</width>
+          <height>0</height>
+         </size>
+        </property>
+        <property name="value">
+         <number>1</number>
+        </property>
+       </widget>
+      </item>
+      <item row="2" column="0">
+       <widget class="QLabel" name="smooth1D_label">
+        <property name="text">
+         <string>Allow 1D relaxation
+ (along borders/constraints)</string>
+        </property>
+        <property name="alignment">
+         <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
+        </property>
+       </widget>
+      </item>
+      <item row="1" column="0">
+       <spacer name="verticalSpacer_3">
+        <property name="orientation">
+         <enum>Qt::Vertical</enum>
+        </property>
+        <property name="sizeType">
+         <enum>QSizePolicy::Maximum</enum>
+        </property>
+        <property name="sizeHint" stdset="0">
+         <size>
+          <width>20</width>
+          <height>40</height>
+         </size>
+        </property>
+       </spacer>
+      </item>
+      <item row="2" column="1">
+       <widget class="QCheckBox" name="smooth1D_checkbox">
+        <property name="text">
+         <string/>
+        </property>
+       </widget>
+      </item>
+     </layout>
+    </widget>
+   </item>
+   <item>
+    <spacer name="verticalSpacer_2">
+     <property name="orientation">
+      <enum>Qt::Vertical</enum>
+     </property>
+     <property name="sizeType">
+      <enum>QSizePolicy::MinimumExpanding</enum>
+     </property>
+     <property name="sizeHint" stdset="0">
+      <size>
+       <width>0</width>
+       <height>0</height>
+      </size>
+     </property>
+    </spacer>
+   </item>
+   <item>
+    <widget class="QDialogButtonBox" name="buttonBox">
+     <property name="orientation">
+      <enum>Qt::Horizontal</enum>
+     </property>
+     <property name="standardButtons">
+      <set>QDialogButtonBox::Cancel|QDialogButtonBox::Ok</set>
+     </property>
+    </widget>
+   </item>
+  </layout>
+ </widget>
+ <tabstops>
+  <tabstop>nbIterations_spinbox</tabstop>
+  <tabstop>smooth1D_checkbox</tabstop>
+ </tabstops>
+ <resources/>
+ <connections>
+  <connection>
+   <sender>buttonBox</sender>
+   <signal>accepted()</signal>
+   <receiver>Tangential_relaxation_dialog</receiver>
+   <slot>accept()</slot>
+   <hints>
+    <hint type="sourcelabel">
+     <x>397</x>
+     <y>333</y>
+    </hint>
+    <hint type="destinationlabel">
+     <x>157</x>
+     <y>195</y>
+    </hint>
+   </hints>
+  </connection>
+  <connection>
+   <sender>buttonBox</sender>
+   <signal>rejected()</signal>
+   <receiver>Tangential_relaxation_dialog</receiver>
+   <slot>reject()</slot>
+   <hints>
+    <hint type="sourcelabel">
+     <x>397</x>
+     <y>333</y>
+    </hint>
+    <hint type="destinationlabel">
+     <x>286</x>
+     <y>195</y>
+    </hint>
+   </hints>
+  </connection>
+ </connections>
+</ui>

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)