From 545009ab7c653107f65eda56251f1f6bde543082 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Sun, 19 Feb 2023 21:25:25 +0000 Subject: [PATCH 1/4] PMP: Add np for not scaling smoothed mesh --- .../internal/Smoothing/curvature_flow_impl.h | 3 ++- .../CGAL/Polygon_mesh_processing/smooth_shape.h | 12 ++++++++++-- .../STL_Extension/internal/parameters_interface.h | 1 + 3 files changed, 13 insertions(+), 3 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h index 411117ab751..0e2978f1d90 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h +++ b/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, + bool scale, const GeomTraits& traits) : mesh_(mesh), vpmap_(vpmap), vcmap_(vcmap), vimap_(get(Vertex_local_index(), mesh_)), - scale_volume_after_smoothing(true), + scale_volume_after_smoothing(scale), traits_(traits), weight_calculator_(mesh_, vpmap_, traits_, false /*no clamping*/, false /*no bounding from below*/) { } diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h index c71cf8b1d0b..d615c67eac1 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h @@ -85,6 +85,13 @@ namespace Polygon_mesh_processing { * \cgalParamExtra{A constrained vertex cannot be modified at all during smoothing.} * \cgalParamNEnd * +* \cgalParamNBegin{scale} +* \cgalParamDescription{A closed mesh with at most one constrained +* vertex can be scaled in order to maintain its volume.} +* \cgalParamType{bool} +* \cgalParamDefault{`true`} +* \cgalParamNEnd +* * \cgalParamNBegin{sparse_linear_solver} * \cgalParamDescription{an instance of the sparse linear solver used for smoothing} * \cgalParamType{a class model of `SparseLinearAlgebraWithFactorTraits_d`} @@ -125,6 +132,7 @@ void smooth_shape(const FaceRange& faces, VCMap vcmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained), Static_boolean_property_map()); const unsigned int nb_iterations = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 1); + const bool scale = choose_parameter(get_parameter(np, internal_np::scale), true); #if defined(CGAL_EIGEN3_ENABLED) #if EIGEN_VERSION_AT_LEAST(3,2,0) @@ -161,11 +169,11 @@ void smooth_shape(const FaceRange& faces, Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n); std::vector > stiffness; - internal::Shape_smoother smoother(tmesh, vpmap, vcmap, gt); + internal::Shape_smoother smoother(tmesh, vpmap, vcmap, scale, 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); diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h index e22b5d6063a..3acbb88fda3 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/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(scale_t, scale, 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) From ee13f778109ca9fc4112843aa2c2b90cda5ddccc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 20 Feb 2023 10:08:21 +0100 Subject: [PATCH 2/4] Various improvements --- .../internal/Smoothing/curvature_flow_impl.h | 16 ++++---- .../Polygon_mesh_processing/smooth_shape.h | 41 +++++++++++-------- .../internal/parameters_interface.h | 2 +- 3 files changed, 32 insertions(+), 27 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h index 0e2978f1d90..9dea9e3b4c1 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h @@ -70,14 +70,14 @@ public: Shape_smoother(TriangleMesh& mesh, VertexPointMap& vpmap, VertexConstraintMap& vcmap, - bool scale, - 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(scale), + scale_volume_after_smoothing_(scale_volume_after_smoothing), traits_(traits), weight_calculator_(mesh_, vpmap_, traits_, false /*no clamping*/, false /*no bounding from below*/) { } @@ -105,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, @@ -223,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; @@ -363,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 anchor_point; // linear system data diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h index d615c67eac1..0d0b9572ca1 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h +++ b/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::%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::%vertex_descriptor` @@ -77,21 +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::%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{scale} -* \cgalParamDescription{A closed mesh with at most one constrained -* vertex can be scaled in order to maintain its volume.} -* \cgalParamType{bool} -* \cgalParamDefault{`true`} -* \cgalParamNEnd -* * \cgalParamNBegin{sparse_linear_solver} * \cgalParamDescription{an instance of the sparse linear solver used for smoothing} * \cgalParamType{a class model of `SparseLinearAlgebraWithFactorTraits_d`} @@ -106,7 +110,8 @@ namespace Polygon_mesh_processing { * * @see `smooth_mesh()` */ -template +template void smooth_shape(const FaceRange& faces, TriangleMesh& tmesh, const double time, @@ -132,7 +137,7 @@ void smooth_shape(const FaceRange& faces, VCMap vcmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained), Static_boolean_property_map()); const unsigned int nb_iterations = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 1); - const bool scale = choose_parameter(get_parameter(np, internal_np::scale), true); + 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) @@ -169,7 +174,7 @@ void smooth_shape(const FaceRange& faces, Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n); std::vector > stiffness; - internal::Shape_smoother smoother(tmesh, vpmap, vcmap, scale, gt); + internal::Shape_smoother smoother(tmesh, vpmap, vcmap, scale_after_smoothing, gt); smoother.init_smoothing(faces); diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h index 3acbb88fda3..2b754467ddf 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -146,7 +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(scale_t, scale, scale) +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) From f2cb368919cecc7a38536231da14f886326c3299 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 20 Feb 2023 10:08:35 +0100 Subject: [PATCH 3/4] Add tests --- .../test_shape_smoothing.cpp | 43 ++++++++++++++++--- 1 file changed, 36 insertions(+), 7 deletions(-) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_shape_smoothing.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_shape_smoothing.cpp index 201d508d085..7ed64301e90 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_shape_smoothing.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_shape_smoothing.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include @@ -17,6 +18,7 @@ #include typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel::FT FT; typedef Kernel::Point_3 Point; typedef CGAL::Surface_mesh SurfaceMesh; typedef CGAL::Polyhedron_3 Mesh_with_id; @@ -38,7 +40,7 @@ void test_implicit_constrained_devil(Mesh mesh) typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typename boost::property_map::type vpmap = get(CGAL::vertex_point, mesh); - // z max is 20 in the devil + // max 'z' is 20 in the devil std::set 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 > vcmap(selected_vertices); std::vector fixed_points(selected_vertices.size()); @@ -124,6 +128,29 @@ void test_implicit_constrained_elephant(Mesh mesh) #endif } +template +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 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(mesh_elephant); test_implicit_constrained_elephant(mesh_elephant); test_implicit_constrained_devil(mesh_devil); + test_implicit_unscaled_elephant(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(pl_mesh_elephant); test_implicit_constrained_elephant(pl_mesh_elephant); test_implicit_constrained_devil(pl_mesh_devil); + test_implicit_unscaled_elephant(pl_mesh_elephant); return EXIT_SUCCESS; } From 4a5ada051a019036b1df562c8258f7d7fcb9f52e Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Thu, 23 Feb 2023 11:35:12 +0000 Subject: [PATCH 4/4] Add to changes.md; Update @see in order to show to non-deprecated function --- Installation/CHANGES.md | 2 ++ .../include/CGAL/Polygon_mesh_processing/smooth_shape.h | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index c5106c6f430..6c255e40af6 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -25,6 +25,8 @@ Release date: June 2023 - 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. diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h index 0d0b9572ca1..bd4e8ad8ab2 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h @@ -108,7 +108,7 @@ 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