Merge pull request #7276 from afabri/PMP_smooth_scale-GF

PMP:  Add np for not scaling smoothed mesh

Installation/CHANGES.md

@@ -34,6 +34,8 @@ CGAL tetrahedral Delaunay refinement algorithm.
 
 -   Added the function `CGAL::Polygon_mesh_processing::remove_almost_degenerate_faces()` to remove badly shaped triangles faces in a mesh.
 
+-   Added a named parameter to `CGAL::Polygon_mesh_processing::smooth_shape()` to disable scaling to compensate volume loss.
+
 ### [3D Simplicial Mesh Data Structure](https://doc.cgal.org/5.6/Manual/packages.html#PkgSMDS3) (new package)
 
 -   This new package wraps all the existing code that deals with a `MeshComplex_3InTriangulation_3` to describe 3D simplicial meshes, and makes the data structure independent from the tetrahedral mesh generation package.

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h

@@ -70,13 +70,14 @@ public:
   Shape_smoother(TriangleMesh& mesh,
                  VertexPointMap& vpmap,
                  VertexConstraintMap& vcmap,
-                 const GeomTraits& traits)
+                 const bool scale_volume_after_smoothing = true,
+                 const GeomTraits& traits = GeomTraits())
     :
       mesh_(mesh),
       vpmap_(vpmap),
       vcmap_(vcmap),
       vimap_(get(Vertex_local_index(), mesh_)),
-      scale_volume_after_smoothing(true),
+      scale_volume_after_smoothing_(scale_volume_after_smoothing),
       traits_(traits),
       weight_calculator_(mesh_, vpmap_, traits_, false /*no clamping*/, false /*no bounding from below*/)
   { }
@@ -104,12 +105,12 @@ public:
         if(anchor_point == boost::none)
           anchor_point = get(vpmap_, v);
         else
-          scale_volume_after_smoothing = false;
+          scale_volume_after_smoothing_ = false;
       }
     }
 
     if(!CGAL::is_closed(mesh_))
-      scale_volume_after_smoothing = false;
+      scale_volume_after_smoothing_ = false;
   }
 
   void setup_system(Eigen_matrix& A,
@@ -222,12 +223,12 @@ public:
   {
     namespace PMP = CGAL::Polygon_mesh_processing;
 
-    if(!scale_volume_after_smoothing)
+    if(!scale_volume_after_smoothing_)
       return update_mesh_no_scaling(Xx, Xy, Xz);
 
     const FT old_vol = volume(mesh_, parameters::vertex_point_map(vpmap_).geom_traits(traits_));
 
-    // If no vertex is constrained, then the smoothed mesh will simply share the same centroid as the input mesh
+    // If no vertex is constrained, then the smoothed mesh will share the same centroid as the input mesh
     Point pre_smooth_anchor_point;
     if(anchor_point != boost::none)
       pre_smooth_anchor_point = *anchor_point;
@@ -362,7 +363,7 @@ private:
   // of volume. We need an anchor point to scale up, either a constrained point or the centroid
   // of the initial mesh if no vertex is constrained. If there is more than a constrained vertex,
   // then no scaling can be done without violating the constraint.
-  bool scale_volume_after_smoothing;
+  bool scale_volume_after_smoothing_;
   boost::optional<Point> anchor_point;
 
   // linear system data

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h

@@ -61,6 +61,25 @@ namespace Polygon_mesh_processing {
 *     \cgalParamDefault{`1`}
 *   \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{do_scale}
+*     \cgalParamDescription{Whether to apply rescaling after smoothing. This is useful because
+*                           the mean curvature flow tends to shrink the surface.}
+*     \cgalParamType{Boolean}
+*     \cgalParamDefault{`true`}
+*     \cgalParamExtra{Scaling can only be applied if the mesh is closed and if there is no more than
+*                     a single constrained vertex.}
+*     \cgalParamExtra{If a vertex is constrained, it is the fixed point of the scaling, otherwise
+*                     the centroid is used.}
+*   \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`
@@ -77,14 +96,6 @@ namespace Polygon_mesh_processing {
 *     \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.}
 *   \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{sparse_linear_solver}
 *     \cgalParamDescription{an instance of the sparse linear solver used for smoothing}
 *     \cgalParamType{a class model of `SparseLinearAlgebraWithFactorTraits_d`}
@@ -97,9 +108,10 @@ namespace Polygon_mesh_processing {
 *
 * @warning This function involves linear algebra, that is computed using non-exact, floating-point arithmetic.
 *
-* @see `smooth_mesh()`
+* @see `angle_and_area_smoothing()`
 */
-template<typename TriangleMesh, typename FaceRange, typename NamedParameters = parameters::Default_named_parameters>
+template<typename TriangleMesh, typename FaceRange,
+         typename NamedParameters = parameters::Default_named_parameters>
 void smooth_shape(const FaceRange& faces,
                   TriangleMesh& tmesh,
                   const double time,
@@ -125,6 +137,7 @@ void smooth_shape(const FaceRange& faces,
   VCMap vcmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained),
                                  Static_boolean_property_map<vertex_descriptor, false>());
   const unsigned int nb_iterations = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 1);
+  const bool scale_after_smoothing = choose_parameter(get_parameter(np, internal_np::do_scale), true);
 
 #if defined(CGAL_EIGEN3_ENABLED)
 #if EIGEN_VERSION_AT_LEAST(3,2,0)
@@ -161,11 +174,11 @@ void smooth_shape(const FaceRange& faces,
   Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n);
   std::vector<CGAL::Triple<std::size_t, std::size_t, double> > stiffness;
 
-  internal::Shape_smoother<TriangleMesh, VertexPointMap, VCMap, Sparse_solver, GeomTraits> smoother(tmesh, vpmap, vcmap, gt);
+  internal::Shape_smoother<TriangleMesh, VertexPointMap, VCMap, Sparse_solver, GeomTraits> smoother(tmesh, vpmap, vcmap, scale_after_smoothing, gt);
 
   smoother.init_smoothing(faces);
 
-  // For robustness reasons, the laplacian coefficients are computed only once (only the mass
+  // For robustness reasons, the Laplacian coefficients are computed only once (only the mass
   // matrix is updated at every iteration). See Kazdhan et al. "Can Mean-Curvature Flow Be Made Non-Singular?".
   smoother.calculate_stiffness_matrix_elements(stiffness);
 

Polygon_mesh_processing/test/Polygon_mesh_processing/test_shape_smoothing.cpp

@@ -6,6 +6,7 @@
 #include <CGAL/Polyhedron_3.h>
 #include <CGAL/Polyhedron_items_with_id_3.h>
 
+#include <CGAL/Polygon_mesh_processing/measure.h>
 #include <CGAL/Polygon_mesh_processing/smooth_shape.h>
 
 #include <CGAL/utility.h>
@@ -17,6 +18,7 @@
 #include <vector>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel                 Kernel;
+typedef Kernel::FT                                                          FT;
 typedef Kernel::Point_3                                                     Point;
 typedef CGAL::Surface_mesh<Point>                                           SurfaceMesh;
 typedef CGAL::Polyhedron_3<Kernel, CGAL::Polyhedron_items_with_id_3>        Mesh_with_id;
@@ -38,7 +40,7 @@ void test_implicit_constrained_devil(Mesh mesh)
   typedef typename boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
   typename boost::property_map<Mesh, CGAL::vertex_point_t>::type vpmap = get(CGAL::vertex_point, mesh);
 
-  // z max is 20 in the devil
+  // max 'z' is 20 in the devil
   std::set<vertex_descriptor> selected_vertices;
   for(vertex_descriptor v : vertices(mesh))
   {
@@ -47,6 +49,8 @@ void test_implicit_constrained_devil(Mesh mesh)
       selected_vertices.insert(v);
   }
 
+  std::cout << selected_vertices.size() << " constrained vertices" << std::endl;
+
   CGAL::Boolean_property_map<std::set<vertex_descriptor> > vcmap(selected_vertices);
 
   std::vector<Point> fixed_points(selected_vertices.size());
@@ -124,6 +128,29 @@ void test_implicit_constrained_elephant(Mesh mesh)
 #endif
 }
 
+template <typename Mesh>
+void test_implicit_unscaled_elephant(Mesh mesh)
+{
+#ifdef CGAL_PMP_SMOOTHING_DEBUG
+  std::cout << "-- test_implicit_unscaled_elephant --" << std::endl;
+#endif
+
+  const FT ivol = PMP::volume(mesh);
+  std::cout << "Input volume is " << ivol << std::endl;
+
+  Mesh mesh_cpy(mesh);
+  const double time_step = 0.001;
+  PMP::smooth_shape(mesh_cpy, time_step, CGAL::parameters::number_of_iterations(5).do_scale(true));
+
+  FT ovol = PMP::volume(mesh_cpy);
+  std::cout << "With scaling, output volume is " << ovol << std::endl;
+  assert(equal_doubles(ivol, ovol, 1e-10));
+
+  PMP::smooth_shape(mesh, time_step, CGAL::parameters::number_of_iterations(5).do_scale(false));
+  ovol = PMP::volume(mesh);
+  std::cout << "Without scaling, output volume is " << ovol << std::endl;
+}
+
 template <typename Mesh>
 void test_curvature_flow_time_step(Mesh mesh)
 {
@@ -167,8 +194,8 @@ int main(int, char**)
   SurfaceMesh mesh_devil;
   if(!input1 || !(input1 >> mesh_devil))
   {
-    std::cerr << "Error: can not read file.";
-    return 1;
+    std::cerr << "Error: cannot read file " << filename_devil << std::endl;
+    return EXIT_FAILURE;
   }
   input1.close();
 
@@ -176,8 +203,8 @@ int main(int, char**)
   SurfaceMesh mesh_elephant;
   if(!input2 || !(input2 >> mesh_elephant))
   {
-    std::cerr << "Error: can not read file.";
-    return 1;
+    std::cerr << "Error: cannot read file " << filename_elephant << std::endl;
+    return EXIT_FAILURE;
   }
   input2.close();
 
@@ -185,11 +212,12 @@ int main(int, char**)
   test_curvature_flow<SurfaceMesh>(mesh_elephant);
   test_implicit_constrained_elephant<SurfaceMesh>(mesh_elephant);
   test_implicit_constrained_devil<SurfaceMesh>(mesh_devil);
+  test_implicit_unscaled_elephant<SurfaceMesh>(mesh_elephant);
 
   input1.open(filename_devil);
   Mesh_with_id pl_mesh_devil;
   if(!input1 || !(input1 >> pl_mesh_devil)){
-    std::cerr << "Error: can not read file.";
+    std::cerr << "Error: cannot read file " << filename_devil << std::endl;
     return EXIT_FAILURE;
   }
   input1.close();
@@ -200,7 +228,7 @@ int main(int, char**)
   Mesh_with_id pl_mesh_elephant;
   if(!input2 || !(input2 >> pl_mesh_elephant))
   {
-    std::cerr << "Error: can not read file.";
+    std::cerr << "Error: cannot read file " << filename_elephant << std::endl;
     return EXIT_FAILURE;
   }
   input2.close();
@@ -212,6 +240,7 @@ int main(int, char**)
   test_curvature_flow<Mesh_with_id>(pl_mesh_elephant);
   test_implicit_constrained_elephant<Mesh_with_id>(pl_mesh_elephant);
   test_implicit_constrained_devil<Mesh_with_id>(pl_mesh_devil);
+  test_implicit_unscaled_elephant<Mesh_with_id>(pl_mesh_elephant);
 
   return EXIT_SUCCESS;
 }

STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h

@@ -146,6 +146,7 @@ CGAL_add_named_parameter(mesh_facet_angle_t, mesh_facet_angle, mesh_facet_angle)
 CGAL_add_named_parameter(mesh_facet_distance_t, mesh_facet_distance, mesh_facet_distance)
 CGAL_add_named_parameter(mesh_facet_topology_t, mesh_facet_topology, mesh_facet_topology)
 CGAL_add_named_parameter(polyline_constraints_t, polyline_constraints, polyline_constraints)
+CGAL_add_named_parameter(do_scale_t, do_scale, do_scale)
 
 // List of named parameters that we use in the package 'Surface Mesh Simplification'
 CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost)