songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

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
This commit is contained in:
Mael Rouxel-Labbé 2019-06-05 15:58:24 +02:00
parent bd7fd4a91e
commit aa6ef78eb3
5 changed files with 136 additions and 133 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_explicit_impl.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2018 GeometryFactory (France). // Copyright (c) 2018-2019 GeometryFactory (France).
// All rights reserved. // All rights reserved.
// //
// This file is part of CGAL (www.cgal.org). // This file is part of CGAL (www.cgal.org).
@ -17,7 +17,8 @@
// SPDX-License-Identifier: GPL-3.0+ // 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 #ifndef CGAL_POLYGON_MESH_PROCESSING_INTERNAL_CURVATURE_FLOW_EXPLICIT_IMPL_H
#define 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 Polygon_mesh_processing {
namespace internal { namespace internal {
// @fixme which cotangent weight should be used?
template<typename TriangleMesh, template<typename TriangleMesh,
typename VertexPointMap, typename VertexPointMap,
typename CotangentValue = CGAL::internal::Cotangent_value_Meyer<TriangleMesh, VertexPointMap> > typename CotangentValue = CGAL::internal::Cotangent_value_Meyer<TriangleMesh, VertexPointMap> >
@ -93,15 +95,18 @@ class Curvature_flow
public: public:
Curvature_flow(TriangleMesh& pmesh, Curvature_flow(TriangleMesh& pmesh,
VertexPointMap vpmap, VertexPointMap vpmap,
VertexConstraintMap vcmap) VertexConstraintMap vcmap,
const GeomTraits& traits)
: :
mesh_(pmesh), mesh_(pmesh),
vpmap_(vpmap), vcmap_(vcmap), vpmap_(vpmap), vcmap_(vcmap),
weight_calculator_(pmesh, vpmap) weight_calculator_(pmesh, vpmap),
traits_(traits)
{ } { }
public:
template<typename FaceRange> template<typename FaceRange>
void init_smoothing(const FaceRange& face_range) void initialize(const FaceRange& face_range)
{ {
CGAL_precondition(CGAL::is_triangle_mesh(mesh_)); CGAL_precondition(CGAL::is_triangle_mesh(mesh_));
CGAL_precondition_code(std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor> degen_faces;) CGAL_precondition_code(std::set<typename boost::graph_traits<TriangleMesh>::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 // @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 // In general, we will smooth whole meshes, so this is better than a map
typedef CGAL::dynamic_vertex_property_t<Point> Vertex_property_tag; typedef CGAL::dynamic_vertex_property_t<Point> Vertex_property_tag;
@ -182,8 +187,6 @@ private:
} }
private: private:
// data members
// ------------
TriangleMesh& mesh_; TriangleMesh& mesh_;
VertexPointMap vpmap_; VertexPointMap vpmap_;
@ -192,7 +195,7 @@ private:
std::vector<vertex_descriptor> vrange_; std::vector<vertex_descriptor> vrange_;
Weight_calculator weight_calculator_; Weight_calculator weight_calculator_;
GeomTraits traits_; const GeomTraits& traits_;
}; };
} // internal } // internal

52
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h
View File

@ -123,6 +123,7 @@ public:
vpmap_(vpmap), vpmap_(vpmap),
vcmap_(vcmap), vcmap_(vcmap),
vimap_(get(Vertex_local_index(), mesh_)), vimap_(get(Vertex_local_index(), mesh_)),
scale_volume_after_smoothing(true),
traits_(traits), traits_(traits),
weight_calculator_(mesh, vpmap) weight_calculator_(mesh, vpmap)
{ } { }
@ -141,18 +142,22 @@ public:
constrained_flags_.assign(vertices(mesh_).size(), false); constrained_flags_.assign(vertices(mesh_).size(), false);
std::size_t counter = 0; // @tmp std::size_t counter = 0; // @tmp
for(vertex_descriptor v : vertices(mesh_))
for(vertex_descriptor v : vrange_)
{ {
if(is_constrained(v)) if(is_constrained(v))
{ {
++counter;
constrained_flags_[get(vimap_, v)] = true; 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, void setup_system(Eigen_matrix& A,
@ -242,10 +247,27 @@ public:
stiffness_elements.push_back(Triplet(it->first, it->first, it->second)); 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) void update_mesh(const Eigen_vector& Xx, const Eigen_vector& Xy, const Eigen_vector& Xz)
{ {
namespace PMP = CGAL::Polygon_mesh_processing; 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_)); 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 simply share the same centroid as the input mesh
@ -314,28 +336,22 @@ private:
{ {
fill_mass_matrix(); fill_mass_matrix();
// fill A = M - time * L // fill A = Mass - time * Laplacian
for(const Triplet& t : stiffness_elements) for(const Triplet& t : stiffness_elements)
{ {
std::size_t i = t.get<0>(), j = t.get<1>(); std::size_t i = t.get<0>(), j = t.get<1>();
if(i != j) if(i != j)
{
A.set_coef(i, j, - time * t.get<2>(), true); 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); A.set_coef(i, i, diagonal_[t.get<0>()] - time * t.get<2>(), true);
}
} }
for(vertex_descriptor v : vrange_) for(vertex_descriptor v : vrange_)
{ {
std::size_t index = get(vimap_, v); std::size_t index = get(vimap_, v);
if(constrained_flags_[index]) if(constrained_flags_[index])
{
std::cout << "A[" << index <<", " << index << "] = 1 (constrained)" << std::endl;
A.set_coef(index, index, 1., true); A.set_coef(index, index, 1., true);
}
} }
// we do not call A.assemble_matrix here // we do not call A.assemble_matrix here
@ -385,9 +401,11 @@ private:
VertexConstraintMap vcmap_; VertexConstraintMap vcmap_;
IndexMap vimap_; IndexMap vimap_;
// Smoothing has a tendency to reduce volumes, so we scale things back up based on the change // Smoothing has a tendency to reduce volumes, so we can 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 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. // 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<Point> anchor_point; boost::optional<Point> anchor_point;
// linear system data // linear system data

3
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h
View File

@ -631,8 +631,7 @@ public:
continue; continue;
Point_ref p_query = get(vpmap_, v); Point_ref p_query = get(vpmap_, v);
Point projected = tree_ptr_->closest_point(p_query); const Point projected = tree_ptr_->closest_point(p_query);
std::cout << p_query << " is projected to: " << projected << std::endl;
put(vpmap_, v, projected); put(vpmap_, v, projected);
} }
} }

34
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_mesh.h
View File

@ -17,8 +17,8 @@
// SPDX-License-Identifier: GPL-3.0+ // SPDX-License-Identifier: GPL-3.0+
// //
// //
// Author(s) : Konstantinos Katrioplas (konst.katrioplas@gmail.com) // Author(s) : Mael Rouxel-Labbé
// Mael Rouxel-Labbé // Konstantinos Katrioplas (konst.katrioplas@gmail.com)
#ifndef CGAL_POLYGON_MESH_PROCESSING_SMOOTH_MESH_H #ifndef CGAL_POLYGON_MESH_PROCESSING_SMOOTH_MESH_H
#define 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<TriangleMesh, VertexPointMap, GeomTraits> Angle_optimizer; typedef internal::Angle_smoother<TriangleMesh, VertexPointMap, GeomTraits> Angle_optimizer;
typedef internal::Mesh_smoother<Angle_optimizer, TriangleMesh, typedef internal::Mesh_smoother<Angle_optimizer, TriangleMesh,
VertexPointMap, VCMap, GeomTraits> Angle_smoother; VertexPointMap, VCMap, GeomTraits> Angle_smoother;
typedef internal::Delaunay_edge_flipper<TriangleMesh, VertexPointMap, GeomTraits> Delaunay_flipper;
if(std::begin(faces) == std::end(faces)) if(std::begin(faces) == std::end(faces))
return; return;
@ -130,10 +129,6 @@ void smooth_angles(const FaceRange& faces,
smoother.project_to_surface(); 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_optimizer, TriangleMesh, typedef internal::Mesh_smoother<Area_optimizer, TriangleMesh,
VertexPointMap, VCMap, GeomTraits> Area_smoother; VertexPointMap, VCMap, GeomTraits> Area_smoother;
typedef internal::Delaunay_edge_flipper<TriangleMesh, VertexPointMap, GeomTraits> Delaunay_flipper;
if(std::begin(faces) == std::end(faces)) if(std::begin(faces) == std::end(faces))
return; return;
@ -248,6 +245,10 @@ void smooth_areas(const FaceRange faces,
smoother.project_to_surface(); 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 <typename FaceRange, typename TriangleMesh> template <typename FaceRange, typename TriangleMesh>
@ -268,15 +269,6 @@ void smooth_areas(TriangleMesh& tmesh)
smooth_areas(faces(tmesh), tmesh, parameters::all_default()); smooth_areas(faces(tmesh), tmesh, parameters::all_default());
} }
///////////////////////////////////////////////////////////////////////////////////////////////////
///
///
///
///
///
///
///////////////////////////////////////////////////////////////////////////////////////////////////
// do both // do both
template<typename TriangleMesh, typename FaceRange, typename NamedParameters> template<typename TriangleMesh, typename FaceRange, typename NamedParameters>
void smooth(const FaceRange& faces, void smooth(const FaceRange& faces,
@ -341,7 +333,11 @@ void smooth(const FaceRange& faces,
area_smoother.project_to_surface(); 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*/, angle_smoother.optimize(use_safety_constraints /*check for bad faces*/,
true /*apply all moves at once*/, true /*apply all moves at once*/,
use_safety_constraints /*check if the min angle is improved*/); use_safety_constraints /*check if the min angle is improved*/);
@ -356,10 +352,6 @@ void smooth(const FaceRange& faces,
angle_smoother.project_to_surface(); 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);
} }
} }

159
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/smooth_shape.h
View File

@ -1,4 +1,4 @@
// Copyright (c) 2018 GeometryFactory (France). // Copyright (c) 2018-2019 GeometryFactory (France).
// All rights reserved. // All rights reserved.
// //
// This file is part of CGAL (www.cgal.org). // This file is part of CGAL (www.cgal.org).
@ -17,7 +17,8 @@
// SPDX-License-Identifier: GPL-3.0+ // 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 #ifndef CGAL_POLYGON_MESH_PROCESSING_SMOOTH_SHAPE_H
#define CGAL_POLYGON_MESH_PROCESSING_SMOOTH_SHAPE_H #define CGAL_POLYGON_MESH_PROCESSING_SMOOTH_SHAPE_H
@ -39,46 +40,35 @@
namespace CGAL { namespace CGAL {
namespace Polygon_mesh_processing { namespace Polygon_mesh_processing {
namespace internal {
// normalized explicit scheme // normalized explicit scheme
template<typename TriangleMesh, typename FaceRange, typename NamedParameters> template<typename FaceRange, typename TriangleMesh,
typename VertexPointMap, typename VertexConstrainMap,
typename GeomTraits>
void smooth_curvature_flow_explicit(const FaceRange& faces, void smooth_curvature_flow_explicit(const FaceRange& faces,
TriangleMesh& tmesh, TriangleMesh& tmesh,
const NamedParameters& np) const VertexPointMap vpmap,
const VertexConstrainMap vcmap,
const std::size_t nb_iterations,
const GeomTraits& traits)
{ {
using boost::choose_param; typedef internal::Curvature_flow<TriangleMesh, VertexPointMap, VertexConstrainMap, GeomTraits> Curvature_explicit_smoother;
using boost::get_param;
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type GeomTraits; Curvature_explicit_smoother curvature_smoother(tmesh, vpmap, vcmap, traits);
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::type VertexPointMap; curvature_smoother.initialize(faces);
VertexPointMap vpmap = choose_param(get_param(np, internal_np::vertex_point),
get_property_map(CGAL::vertex_point, tmesh));
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::lookup_named_param_def <
internal_np::vertex_is_constrained_t,
NamedParameters,
Constant_property_map<vertex_descriptor, bool>
> ::type VCMap;
VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained),
Constant_property_map<vertex_descriptor, bool>(false));
std::size_t nb_iterations = choose_param(get_param(np, internal_np::number_of_iterations), 1);
internal::Curvature_flow<TriangleMesh, VertexPointMap, VCMap, GeomTraits>
curvature_smoother(tmesh, vpmap, vcmap);
curvature_smoother.init_smoothing(faces);
for(std::size_t i=0; i<nb_iterations; ++i) for(std::size_t i=0; i<nb_iterations; ++i)
{ {
#ifdef CGAL_PMP_SMOOTHING_VERBOSE #ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cout << "iteration #" << i << std::endl; std::cout << "iteration #" << i << std::endl;
#endif #endif
curvature_smoother.curvature_smoothing(); curvature_smoother.smooth();
} }
} }
} // namespace internal
/*! /*!
* \ingroup PMP_meshing_grp * \ingroup PMP_meshing_grp
* smooths the overall shape of the mesh by using the mean curvature flow. * smooths the overall shape of the mesh by using the mean curvature flow.
@ -132,94 +122,95 @@ void smooth_along_curvature_flow(const FaceRange& faces,
if(std::begin(faces) == std::end(faces)) if(std::begin(faces) == std::end(faces))
return; return;
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type GeomTraits;
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::type VertexPointMap;
typedef typename boost::lookup_named_param_def<
internal_np::vertex_is_constrained_t,
NamedParameters,
Constant_property_map<vertex_descriptor, bool> >::type VCMap;
using boost::choose_param; using boost::choose_param;
using boost::get_param; using boost::get_param;
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type GeomTraits;
GeomTraits gt = choose_param(get_param(np, internal_np::geom_traits), GeomTraits()); GeomTraits gt = choose_param(get_param(np, internal_np::geom_traits), GeomTraits());
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::type VertexPointMap;
VertexPointMap vpmap = choose_param(get_param(np, internal_np::vertex_point), VertexPointMap vpmap = choose_param(get_param(np, internal_np::vertex_point),
get_property_map(CGAL::vertex_point, tmesh)); get_property_map(CGAL::vertex_point, tmesh));
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::lookup_named_param_def<
internal_np::vertex_is_constrained_t,
NamedParameters,
Constant_property_map<vertex_descriptor, bool>
> ::type VCMap;
VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained), VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained),
Constant_property_map<vertex_descriptor, bool>(false)); Constant_property_map<vertex_descriptor, bool>(false));
const std::size_t nb_iterations = choose_param(get_param(np, internal_np::number_of_iterations), 1);
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);
bool use_explicit_scheme = choose_param(get_param(np, internal_np::use_explicit_scheme), false);
if(use_explicit_scheme) if(use_explicit_scheme)
{ return internal::smooth_curvature_flow_explicit(faces, tmesh, vpmap, vcmap, nb_iterations, gt);
smooth_curvature_flow_explicit(faces, tmesh, parameters::number_of_iterations(nb_iterations));
} // from now on, it is the implicit scheme
else
{
// implicit scheme
#if defined(CGAL_EIGEN3_ENABLED) #if defined(CGAL_EIGEN3_ENABLED)
#if EIGEN_VERSION_AT_LEAST(3,2,0) #if EIGEN_VERSION_AT_LEAST(3,2,0)
typedef typename Eigen::SparseMatrix<double> Eigen_sparse_matrix; typedef typename Eigen::SparseMatrix<double> Eigen_sparse_matrix;
typedef typename Eigen::BiCGSTAB<Eigen_sparse_matrix, Eigen::IncompleteLUT<double> > Eigen_solver; typedef typename Eigen::BiCGSTAB<Eigen_sparse_matrix, Eigen::IncompleteLUT<double> > Eigen_solver;
typedef CGAL::Eigen_solver_traits<Eigen_solver> Default_solver; typedef CGAL::Eigen_solver_traits<Eigen_solver> Default_solver;
#else
typedef bool Default_solver;//compilation should crash
//if no solver is provided and Eigen version < 3.2
#endif
#else #else
typedef bool Default_solver;//compilation should crash typedef bool Default_solver;//compilation should crash
// if no solver is provided and Eigen version < 3.2 //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 #endif
#if defined(CGAL_EIGEN3_ENABLED) #if defined(CGAL_EIGEN3_ENABLED)
CGAL_static_assertion_msg( CGAL_static_assertion_msg(
(!boost::is_same<typename GetSolver<NamedParameters, Default_solver>::type, bool>::value) || EIGEN_VERSION_AT_LEAST(3, 2, 0), (!boost::is_same<typename GetSolver<NamedParameters, Default_solver>::type, bool>::value) || EIGEN_VERSION_AT_LEAST(3, 2, 0),
"Eigen3 version 3.2 or later is required."); "Eigen3 version 3.2 or later is required.");
#else #else
CGAL_static_assertion_msg( CGAL_static_assertion_msg(
(!boost::is_same<typename GetSolver<NamedParameters, Default_solver>::type, bool>::value), (!boost::is_same<typename GetSolver<NamedParameters, Default_solver>::type, bool>::value),
"Eigen3 version 3.2 or later is required."); "Eigen3 version 3.2 or later is required.");
#endif #endif
typedef typename GetSolver<NamedParameters, Default_solver>::type Sparse_solver; typedef typename GetSolver<NamedParameters, Default_solver>::type Sparse_solver;
typedef typename Sparse_solver::Matrix Eigen_matrix; typedef typename Sparse_solver::Matrix Eigen_matrix;
typedef typename Sparse_solver::Vector Eigen_vector; 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(); std::size_t n = vertices(tmesh).size();
Eigen_matrix A(n, n); Eigen_matrix A(n, n);
Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n); 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; 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, gt);
smoother.init_smoothing(faces); 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?". // matrix is updated at every iteration). See Kazdhan et al. "Can Mean-Curvature Flow Be Made Non-Singular?".
smoother.calculate_stiffness_matrix_elements(stiffness); smoother.calculate_stiffness_matrix_elements(stiffness);
for(std::size_t iter = 0; iter < nb_iterations; ++iter) for(std::size_t iter=0; iter<nb_iterations; ++iter)
{ {
std::cout << "iter: " << iter << std::endl; #ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::cout << "iteration #" << i << std::endl;
#endif
smoother.setup_system(A, bx, by, bz, stiffness, time); smoother.setup_system(A, bx, by, bz, stiffness, time);
smoother.solve_system(A, Xx, Xy, Xz, bx, by, bz, solver);
if(smoother.solve_system(A, Xx, Xy, Xz, bx, by, bz, solver))
smoother.update_mesh(Xx, Xy, Xz); smoother.update_mesh(Xx, Xy, Xz);
else
break;
// @tmp (clean fstream sstream includes too) #ifdef CGAL_PMP_SMOOTHING_VERBOSE
std::stringstream oss; // @tmp (clean fstream sstream includes too)
oss << "intermediate_" << iter << ".off" << std::ends; std::stringstream oss;
std::ofstream out(oss.str().c_str()); oss << "intermediate_" << iter << ".off" << std::ends;
out.precision(17); std::ofstream out(oss.str().c_str());
out << tmesh; out.precision(17);
out.close(); out << tmesh;
} out.close();
#endif
} }
} }