From aa6ef78eb35ff0f5bb22cedb78d2d105b6713ef1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Wed, 5 Jun 2019 15:58:24 +0200 Subject: [PATCH] Various fixes/improvements: - Fix using Delaunay flips after angle smoothing (must be after area smoothing) - Fix parameter passing to explicit curvature flow formulation - Use early exits - Clarify function names --- .../Smoothing/curvature_flow_explicit_impl.h | 21 ++- .../internal/Smoothing/curvature_flow_impl.h | 52 ++++-- .../internal/Smoothing/mesh_smoothing_impl.h | 3 +- .../Polygon_mesh_processing/smooth_mesh.h | 34 ++-- .../Polygon_mesh_processing/smooth_shape.h | 159 +++++++++--------- 5 files changed, 136 insertions(+), 133 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_explicit_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_explicit_impl.h index b879124bcab..88cf6023c51 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_explicit_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_explicit_impl.h @@ -1,4 +1,4 @@ -// Copyright (c) 2018 GeometryFactory (France). +// Copyright (c) 2018-2019 GeometryFactory (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -17,7 +17,8 @@ // SPDX-License-Identifier: GPL-3.0+ // // -// Author(s) : Konstantinos Katrioplas (konst.katrioplas@gmail.com) +// Author(s) : Mael Rouxel-Labbé +// Konstantinos Katrioplas (konst.katrioplas@gmail.com) #ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CURVATURE_FLOW_EXPLICIT_IMPL_H #define CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CURVATURE_FLOW_EXPLICIT_IMPL_H @@ -46,6 +47,7 @@ namespace CGAL { namespace Polygon_mesh_processing { namespace internal { +// @fixme which cotangent weight should be used? template > @@ -93,15 +95,18 @@ class Curvature_flow public: Curvature_flow(TriangleMesh& pmesh, VertexPointMap vpmap, - VertexConstraintMap vcmap) + VertexConstraintMap vcmap, + const GeomTraits& traits) : mesh_(pmesh), vpmap_(vpmap), vcmap_(vcmap), - weight_calculator_(pmesh, vpmap) + weight_calculator_(pmesh, vpmap), + traits_(traits) { } +public: template - void init_smoothing(const FaceRange& face_range) + void initialize(const FaceRange& face_range) { CGAL_precondition(CGAL::is_triangle_mesh(mesh_)); CGAL_precondition_code(std::set::face_descriptor> degen_faces;) @@ -114,7 +119,7 @@ public: } // @todo should be possible to pass a time step because this is not a very stable process - void curvature_smoothing() + void smooth() { // In general, we will smooth whole meshes, so this is better than a map typedef CGAL::dynamic_vertex_property_t Vertex_property_tag; @@ -182,8 +187,6 @@ private: } private: - // data members - // ------------ TriangleMesh& mesh_; VertexPointMap vpmap_; @@ -192,7 +195,7 @@ private: std::vector vrange_; Weight_calculator weight_calculator_; - GeomTraits traits_; + const GeomTraits& traits_; }; } // internal 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 132b4b66e3c..bad96c2d3b1 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 @@ -123,6 +123,7 @@ public: vpmap_(vpmap), vcmap_(vcmap), vimap_(get(Vertex_local_index(), mesh_)), + scale_volume_after_smoothing(true), traits_(traits), weight_calculator_(mesh, vpmap) { } @@ -141,18 +142,22 @@ public: constrained_flags_.assign(vertices(mesh_).size(), false); std::size_t counter = 0; // @tmp - - for(vertex_descriptor v : vrange_) + for(vertex_descriptor v : vertices(mesh_)) { if(is_constrained(v)) { - ++counter; constrained_flags_[get(vimap_, v)] = true; - anchor_point = get(vpmap_, v); + ++counter; + + // scaling things cannot preserve the position of more than a single constrained point + if(anchor_point == boost::none) + anchor_point = get(vpmap_, v); + else + scale_volume_after_smoothing = false; } } - std::cout << counter << " vertices are constrained" << std::endl; + std::cout << counter << " constrained vertices" << std::endl; } void setup_system(Eigen_matrix& A, @@ -242,10 +247,27 @@ public: stiffness_elements.push_back(Triplet(it->first, it->first, it->second)); } + void update_mesh_no_scaling(const Eigen_vector& Xx, const Eigen_vector& Xy, const Eigen_vector& Xz) + { + for(vertex_descriptor v : vrange_) + { + std::size_t index = get(vimap_, v); + const FT x_new = Xx[index]; + const FT y_new = Xy[index]; + const FT z_new = Xz[index]; + + Point new_pos(x_new, y_new, z_new); + put(vpmap_, v, new_pos); + } + } + void update_mesh(const Eigen_vector& Xx, const Eigen_vector& Xy, const Eigen_vector& Xz) { namespace PMP = CGAL::Polygon_mesh_processing; + 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 @@ -314,28 +336,22 @@ private: { fill_mass_matrix(); - // fill A = M - time * L + // fill A = Mass - time * Laplacian for(const Triplet& t : stiffness_elements) { std::size_t i = t.get<0>(), j = t.get<1>(); + if(i != j) - { A.set_coef(i, j, - time * t.get<2>(), true); - } - else if(!constrained_flags_[i]) // i==j - { + else if(!constrained_flags_[i]) // && i==j A.set_coef(i, i, diagonal_[t.get<0>()] - time * t.get<2>(), true); - } } for(vertex_descriptor v : vrange_) { std::size_t index = get(vimap_, v); if(constrained_flags_[index]) - { - std::cout << "A[" << index <<", " << index << "] = 1 (constrained)" << std::endl; A.set_coef(index, index, 1., true); - } } // we do not call A.assemble_matrix here @@ -385,9 +401,11 @@ private: VertexConstraintMap vcmap_; IndexMap vimap_; - // Smoothing has a tendency to reduce volumes, so we scale things back up based on the change - // of volume. We also need an anchor point to scale up, either a constrained point or the centroid - // of the initial mesh if no vertex is constrained. + // Smoothing has a tendency to reduce volumes, so we can scale things back up based on the change + // 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; boost::optional anchor_point; // linear system data diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h index bd106e054d7..c3253d74ee2 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h @@ -631,8 +631,7 @@ public: continue; Point_ref p_query = get(vpmap_, v); - Point projected = tree_ptr_->closest_point(p_query); - std::cout << p_query << " is projected to: " << projected << std::endl; + const Point projected = tree_ptr_->closest_point(p_query); put(vpmap_, v, projected); } } diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h index 8ee12112f40..407b2dc087d 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h @@ -17,8 +17,8 @@ // SPDX-License-Identifier: GPL-3.0+ // // -// Author(s) : Konstantinos Katrioplas (konst.katrioplas@gmail.com) -// Mael Rouxel-Labbé +// Author(s) : Mael Rouxel-Labbé +// Konstantinos Katrioplas (konst.katrioplas@gmail.com) #ifndef CGAL_POLYGON_MESH_PROCESSING_SMOOTH_MESH_H #define CGAL_POLYGON_MESH_PROCESSING_SMOOTH_MESH_H @@ -90,7 +90,6 @@ void smooth_angles(const FaceRange& faces, typedef internal::Angle_smoother Angle_optimizer; typedef internal::Mesh_smoother Angle_smoother; - typedef internal::Delaunay_edge_flipper Delaunay_flipper; if(std::begin(faces) == std::end(faces)) return; @@ -130,10 +129,6 @@ void smooth_angles(const FaceRange& faces, smoother.project_to_surface(); } - - // according to the paper, we're supposed to Delaunay-based edge flips! - Delaunay_flipper delaunay_flipper(tmesh, vpmap, gt); - delaunay_flipper(faces); } } @@ -209,6 +204,8 @@ void smooth_areas(const FaceRange faces, typedef internal::Mesh_smoother Area_smoother; + typedef internal::Delaunay_edge_flipper Delaunay_flipper; + if(std::begin(faces) == std::end(faces)) return; @@ -248,6 +245,10 @@ void smooth_areas(const FaceRange faces, smoother.project_to_surface(); } } + + // according to the paper, we're supposed to Delaunay-based edge flips! + Delaunay_flipper delaunay_flipper(tmesh, vpmap, gt); + delaunay_flipper(faces); } template @@ -268,15 +269,6 @@ void smooth_areas(TriangleMesh& tmesh) smooth_areas(faces(tmesh), tmesh, parameters::all_default()); } -/////////////////////////////////////////////////////////////////////////////////////////////////// -/// -/// -/// -/// -/// -/// -/////////////////////////////////////////////////////////////////////////////////////////////////// - // do both template void smooth(const FaceRange& faces, @@ -341,7 +333,11 @@ void smooth(const FaceRange& faces, area_smoother.project_to_surface(); } - // ... then angle smoothing + Delaunay flip + // @todo this should be optional + Delaunay_flipper delaunay_flipper(tmesh, vpmap, gt); + delaunay_flipper(faces); + + // ... then angle smoothing 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*/); @@ -356,10 +352,6 @@ void smooth(const FaceRange& faces, angle_smoother.project_to_surface(); } - - // according to the paper, we're supposed to Delaunay-based edge flips! - Delaunay_flipper delaunay_flipper(tmesh, vpmap, gt); - delaunay_flipper(faces); } } 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 a6935440483..53e983cc449 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 @@ -1,4 +1,4 @@ -// Copyright (c) 2018 GeometryFactory (France). +// Copyright (c) 2018-2019 GeometryFactory (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -17,7 +17,8 @@ // SPDX-License-Identifier: GPL-3.0+ // // -// Author(s) : Konstantinos Katrioplas (konst.katrioplas@gmail.com) +// Author(s) : Mael Rouxel-Labbé +// Konstantinos Katrioplas (konst.katrioplas@gmail.com) #ifndef CGAL_POLYGON_MESH_PROCESSING_SMOOTH_SHAPE_H #define CGAL_POLYGON_MESH_PROCESSING_SMOOTH_SHAPE_H @@ -39,46 +40,35 @@ namespace CGAL { namespace Polygon_mesh_processing { +namespace internal { // normalized explicit scheme -template +template void smooth_curvature_flow_explicit(const FaceRange& faces, TriangleMesh& tmesh, - const NamedParameters& np) + const VertexPointMap vpmap, + const VertexConstrainMap vcmap, + const std::size_t nb_iterations, + const GeomTraits& traits) { - using boost::choose_param; - using boost::get_param; + typedef internal::Curvature_flow Curvature_explicit_smoother; - typedef typename GetGeomTraits::type GeomTraits; - typedef typename GetVertexPointMap::type VertexPointMap; - VertexPointMap vpmap = choose_param(get_param(np, internal_np::vertex_point), - get_property_map(CGAL::vertex_point, tmesh)); - - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::lookup_named_param_def < - internal_np::vertex_is_constrained_t, - NamedParameters, - Constant_property_map - > ::type VCMap; - VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained), - Constant_property_map(false)); - - std::size_t nb_iterations = choose_param(get_param(np, internal_np::number_of_iterations), 1); - - internal::Curvature_flow - curvature_smoother(tmesh, vpmap, vcmap); - - curvature_smoother.init_smoothing(faces); + Curvature_explicit_smoother curvature_smoother(tmesh, vpmap, vcmap, traits); + curvature_smoother.initialize(faces); for(std::size_t i=0; i::vertex_descriptor vertex_descriptor; + typedef typename GetGeomTraits::type GeomTraits; + typedef typename GetVertexPointMap::type VertexPointMap; + typedef typename boost::lookup_named_param_def< + internal_np::vertex_is_constrained_t, + NamedParameters, + Constant_property_map >::type VCMap; + using boost::choose_param; using boost::get_param; - typedef typename GetGeomTraits::type GeomTraits; GeomTraits gt = choose_param(get_param(np, internal_np::geom_traits), GeomTraits()); - - typedef typename GetVertexPointMap::type VertexPointMap; VertexPointMap vpmap = choose_param(get_param(np, internal_np::vertex_point), get_property_map(CGAL::vertex_point, tmesh)); - - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::lookup_named_param_def< - internal_np::vertex_is_constrained_t, - NamedParameters, - Constant_property_map - > ::type VCMap; VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained), Constant_property_map(false)); - - std::size_t nb_iterations = choose_param(get_param(np, internal_np::number_of_iterations), 1); - bool use_explicit_scheme = choose_param(get_param(np, internal_np::use_explicit_scheme), false); + const std::size_t nb_iterations = choose_param(get_param(np, internal_np::number_of_iterations), 1); + const bool use_explicit_scheme = choose_param(get_param(np, internal_np::use_explicit_scheme), false); if(use_explicit_scheme) - { - smooth_curvature_flow_explicit(faces, tmesh, parameters::number_of_iterations(nb_iterations)); - } - else - { - // implicit scheme + return internal::smooth_curvature_flow_explicit(faces, tmesh, vpmap, vcmap, nb_iterations, gt); + + // from now on, it is the implicit scheme + #if defined(CGAL_EIGEN3_ENABLED) - #if EIGEN_VERSION_AT_LEAST(3,2,0) - typedef typename Eigen::SparseMatrix Eigen_sparse_matrix; - typedef typename Eigen::BiCGSTAB > Eigen_solver; - typedef CGAL::Eigen_solver_traits Default_solver; - #else - typedef bool Default_solver;//compilation should crash - //if no solver is provided and Eigen version < 3.2 - #endif +#if EIGEN_VERSION_AT_LEAST(3,2,0) + typedef typename Eigen::SparseMatrix Eigen_sparse_matrix; + typedef typename Eigen::BiCGSTAB > Eigen_solver; + typedef CGAL::Eigen_solver_traits Default_solver; #else - typedef bool Default_solver;//compilation should crash - // if no solver is provided and Eigen version < 3.2 + typedef bool Default_solver;//compilation should crash + //if no solver is provided and Eigen version < 3.2 +#endif +#else + typedef bool Default_solver;//compilation should crash + // if no solver is provided and Eigen version < 3.2 #endif #if defined(CGAL_EIGEN3_ENABLED) - CGAL_static_assertion_msg( + CGAL_static_assertion_msg( (!boost::is_same::type, bool>::value) || EIGEN_VERSION_AT_LEAST(3, 2, 0), "Eigen3 version 3.2 or later is required."); #else - CGAL_static_assertion_msg( + CGAL_static_assertion_msg( (!boost::is_same::type, bool>::value), "Eigen3 version 3.2 or later is required."); #endif - typedef typename GetSolver::type Sparse_solver; - typedef typename Sparse_solver::Matrix Eigen_matrix; - typedef typename Sparse_solver::Vector Eigen_vector; + typedef typename GetSolver::type Sparse_solver; + typedef typename Sparse_solver::Matrix Eigen_matrix; + typedef typename Sparse_solver::Vector Eigen_vector; - Sparse_solver solver = choose_param(get_param(np, internal_np::sparse_linear_solver), Default_solver()); + Sparse_solver solver = choose_param(get_param(np, internal_np::sparse_linear_solver), Default_solver()); - std::size_t n = vertices(tmesh).size(); - Eigen_matrix A(n, n); - Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n); - std::vector > stiffness; + std::size_t n = vertices(tmesh).size(); + Eigen_matrix A(n, n); + 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, gt); - smoother.init_smoothing(faces); + smoother.init_smoothing(faces); - // 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); + // 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); - for(std::size_t iter = 0; iter < nb_iterations; ++iter) - { - std::cout << "iter: " << iter << std::endl; + for(std::size_t iter=0; iter