mirror of https://github.com/CGAL/cgal
add implicit scheme for shape smoothing
This commit is contained in:
konstantinos katrioplas
parent
52126819e0
commit
23670bb712
|
|
@ -0,0 +1,627 @@
|
|||
#ifndef SMOOTHING_H
|
||||
#define SMOOTHING_H
|
||||
|
||||
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
|
||||
#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
|
||||
#include <CGAL/Polygon_mesh_processing/Weights.h>
|
||||
#include <CGAL/Polygon_mesh_processing/measure.h>
|
||||
#include <boost/graph/graph_traits.hpp>
|
||||
#include <boost/property_map/property_map.hpp>
|
||||
|
||||
#include <CGAL/Eigen_solver_traits.h>
|
||||
#include <CGAL/Eigen_matrix.h>
|
||||
#include <CGAL/barycenter.h>
|
||||
#include <Eigen/Sparse>
|
||||
|
||||
#include <fstream>
|
||||
#include <unordered_set>
|
||||
#include <unordered_map>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
namespace Polygon_mesh_processing {
|
||||
|
||||
namespace internal {
|
||||
//Cotangent_value_Meyer<PolygonMesh, VertexPointMap>>
|
||||
|
||||
std::ostream& nl(std::ostream& out)
|
||||
{
|
||||
return out << "\n";
|
||||
}
|
||||
|
||||
template<typename PolygonMesh, typename VertexPointMap,
|
||||
typename CotangentValue = CGAL::internal::Cotangent_value_Meyer<PolygonMesh, VertexPointMap>>
|
||||
struct Edge_cotangent_weight : CotangentValue
|
||||
{
|
||||
Edge_cotangent_weight(PolygonMesh& pmesh_, VertexPointMap vpmap_)
|
||||
: CotangentValue(pmesh_, vpmap_)
|
||||
{}
|
||||
|
||||
PolygonMesh& pmesh()
|
||||
{
|
||||
return CotangentValue::pmesh();
|
||||
}
|
||||
|
||||
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
|
||||
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
|
||||
|
||||
double operator()(halfedge_descriptor he)
|
||||
{
|
||||
if(is_border_edge(he, pmesh()))
|
||||
{
|
||||
halfedge_descriptor h1 = next(he, pmesh());
|
||||
vertex_descriptor vs = source(he, pmesh());
|
||||
vertex_descriptor vt = target(he, pmesh());
|
||||
vertex_descriptor v1 = target(h1, pmesh());
|
||||
|
||||
return (CotangentValue::operator ()(vs, v1, vt));
|
||||
}
|
||||
else
|
||||
{
|
||||
halfedge_descriptor h1 = next(he, pmesh());
|
||||
halfedge_descriptor h2 = prev(opposite(he, pmesh()), pmesh());
|
||||
vertex_descriptor vs = source(he, pmesh());
|
||||
vertex_descriptor vt = target(he, pmesh());
|
||||
vertex_descriptor v1 = target(h1, pmesh());
|
||||
vertex_descriptor v2 = source(h2, pmesh());
|
||||
|
||||
return ( CotangentValue::operator()(vs, v1, vt) + CotangentValue::operator()(vs, v2, vt) ) / 2.0;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename PolygonMesh>
|
||||
struct Incident_area
|
||||
{
|
||||
|
||||
Incident_area(PolygonMesh& mesh) : pmesh(mesh){}
|
||||
|
||||
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
|
||||
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
|
||||
|
||||
double operator()(halfedge_descriptor he)
|
||||
{
|
||||
|
||||
halfedge_descriptor hopp = opposite(he, pmesh);
|
||||
face_descriptor f1 = face(he, pmesh);
|
||||
face_descriptor f2 = face(hopp, pmesh);
|
||||
|
||||
double A1 = f1 == boost::graph_traits<PolygonMesh>::null_face() ? 0 : face_area(f1, pmesh);
|
||||
double A2 = f2 == boost::graph_traits<PolygonMesh>::null_face() ? 0 : face_area(f2, pmesh);
|
||||
|
||||
// todo: check degenerecies
|
||||
//CGAL_assertion(A1>0 && A2>0);
|
||||
|
||||
return A1 + A2;
|
||||
}
|
||||
|
||||
PolygonMesh& pmesh;
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename PolygonMesh, typename VertexPointMap>
|
||||
class Shape_smoother{
|
||||
|
||||
// types
|
||||
private:
|
||||
|
||||
typedef typename GetGeomTraits<PolygonMesh>::type GeomTraits;
|
||||
typedef typename GeomTraits::FT NT;
|
||||
typedef typename GeomTraits::Point_3 Point;
|
||||
typedef typename GeomTraits::Triangle_3 Triangle;
|
||||
|
||||
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
|
||||
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
|
||||
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
|
||||
|
||||
// vertex index map
|
||||
typedef typename boost::property_map<PolygonMesh, boost::vertex_index_t>::type IndexMap;
|
||||
|
||||
//CGAL Sparse solver
|
||||
typedef typename CGAL::Eigen_sparse_matrix<double>::EigenType CEigenMatrix;
|
||||
|
||||
typedef CGAL::Eigen_solver_traits< Eigen::SimplicialLDLT< CEigenMatrix > > Solver_traits;
|
||||
//typedef CGAL::Eigen_solver_traits<> Solver_traits; //BicGSTAB
|
||||
|
||||
typedef typename Solver_traits::Matrix Matrix;
|
||||
typedef typename Solver_traits::Vector Vector;
|
||||
|
||||
// Eigen sparse matrix & vector
|
||||
typedef typename Eigen::SparseMatrix<double> Eigen_matrix;
|
||||
typedef typename Eigen::SparseVector<double> Eigen_vector;
|
||||
|
||||
|
||||
|
||||
// data
|
||||
private:
|
||||
|
||||
IndexMap vimap_ = get(boost::vertex_index, mesh_);
|
||||
|
||||
// geometry data
|
||||
PolygonMesh& mesh_;
|
||||
VertexPointMap& vpmap_;
|
||||
Edge_cotangent_weight<PolygonMesh, VertexPointMap> weight_calculator_;
|
||||
Incident_area<PolygonMesh> inc_areas_calculator_;
|
||||
|
||||
std::size_t nb_vert_;
|
||||
|
||||
// linear solver
|
||||
Solver_traits solver_;
|
||||
|
||||
// constrained vertices
|
||||
std::unordered_set<vertex_descriptor> constrained_vertices_;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
Shape_smoother(PolygonMesh& mesh, VertexPointMap& vpmap) : mesh_(mesh), vpmap_(vpmap),
|
||||
weight_calculator_(mesh, vpmap),
|
||||
inc_areas_calculator_(mesh),
|
||||
nb_vert_(static_cast<int>(vertices(mesh).size()))
|
||||
{ }
|
||||
|
||||
|
||||
Eigen_matrix calc_stiff_matrix()
|
||||
{
|
||||
Eigen_matrix L = stiff_matrix();
|
||||
return L;
|
||||
}
|
||||
|
||||
|
||||
void solve_system(Eigen::SparseMatrix<double>& L)
|
||||
{
|
||||
|
||||
//gather_constrained_vertices();
|
||||
|
||||
Matrix A(nb_vert_, nb_vert_);
|
||||
Vector Bx(nb_vert_);
|
||||
Vector By(nb_vert_);
|
||||
Vector Bz(nb_vert_);
|
||||
Vector Xx(nb_vert_);
|
||||
Vector Xy(nb_vert_);
|
||||
Vector Xz(nb_vert_);
|
||||
|
||||
|
||||
compute_coeff_matrix(A, L);
|
||||
|
||||
//apply_constraints(A);
|
||||
|
||||
|
||||
//extract_matrix(A);
|
||||
|
||||
|
||||
compute_rhs(Bx, By, Bz);
|
||||
|
||||
//extract_vectors(Bx,By,Bz);
|
||||
|
||||
|
||||
NT dx, dy, dz;
|
||||
if(!solver_.linear_solver(A, Bx, Xx, dx) ||
|
||||
!solver_.linear_solver(A, By, Xy, dy) ||
|
||||
!solver_.linear_solver(A, Bz, Xz, dz) )
|
||||
{
|
||||
std::cerr<<"Could not solve linear system!"<<nl;
|
||||
bool solved = false;
|
||||
CGAL_assertion(solved);
|
||||
}
|
||||
|
||||
//extract_vectors(Xx,Xy,Xz);
|
||||
|
||||
//translate_centroid(Xx, Xy, Xz);
|
||||
//extract_vectors(Xx,Xy,Xz);
|
||||
|
||||
NT surface_area = area(faces(mesh_), mesh_);
|
||||
std::cout << "area= " << surface_area << nl;
|
||||
std::cout << "area sqrt= " << CGAL::sqrt(surface_area) << nl;
|
||||
|
||||
|
||||
normalize_area(Xx, Xy, Xz);
|
||||
//extract_vectors(Xx, Xy, Xz);
|
||||
|
||||
update_map(Xx, Xy, Xz);
|
||||
|
||||
surface_area = area(faces(mesh_), mesh_);
|
||||
std::cout << "surface_area normalized= " << surface_area << nl;
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
private:
|
||||
|
||||
// ----------- LINEAR SYSTEM -------- //
|
||||
Eigen_matrix get_stiffness_matrix()
|
||||
{
|
||||
Eigen_matrix mat(nb_vert_, nb_vert_);
|
||||
for(vertex_descriptor vi : vertices(mesh_))
|
||||
{
|
||||
NT sum_Lik = 0;
|
||||
for(halfedge_descriptor h : halfedges_around_source(vi, mesh_))
|
||||
{
|
||||
NT Lij = weight_calculator_(h);
|
||||
sum_Lik -= Lij;
|
||||
vertex_descriptor vj = target(h, mesh_);
|
||||
mat.coeffRef(vimap_[vi], vimap_[vj]) = Lij;
|
||||
}
|
||||
|
||||
mat.coeffRef(vimap_[vi], vimap_[vi]) = sum_Lik;
|
||||
}
|
||||
|
||||
return mat;
|
||||
}
|
||||
|
||||
Eigen_matrix get_mass_matrix()
|
||||
{
|
||||
Eigen_matrix mat(nb_vert_, nb_vert_);
|
||||
for(vertex_descriptor vi : vertices(mesh_))
|
||||
{
|
||||
NT sum_Dik = 0;
|
||||
for(halfedge_descriptor h : halfedges_around_source(vi, mesh_))
|
||||
{
|
||||
NT Dij = inc_areas_calculator_(h) / 12.0;
|
||||
sum_Dik += Dij;
|
||||
vertex_descriptor vj = target(h, mesh_);
|
||||
mat.coeffRef(vimap_[vi], vimap_[vj]) = Dij;
|
||||
}
|
||||
|
||||
mat.coeffRef(vimap_[vi], vimap_[vi]) = sum_Dik;
|
||||
}
|
||||
|
||||
return mat;
|
||||
}
|
||||
|
||||
|
||||
Eigen_matrix cot_entries()
|
||||
{
|
||||
Eigen_matrix C(nb_vert_, nb_vert_);
|
||||
|
||||
for(vertex_descriptor v : vertices(mesh_))
|
||||
{
|
||||
for(halfedge_descriptor h : halfedges_around_source(v, mesh_))
|
||||
{
|
||||
NT Lij = weight_calculator_(h);
|
||||
|
||||
vertex_descriptor v_source = source(h, mesh_);
|
||||
vertex_descriptor v_target = target(h, mesh_);
|
||||
auto i_source = vimap_[v_source];
|
||||
auto i_target = vimap_[v_target];
|
||||
|
||||
C.coeffRef(i_source, i_target) = Lij;
|
||||
}
|
||||
}
|
||||
|
||||
return C;
|
||||
}
|
||||
|
||||
Eigen_matrix stiff_matrix()
|
||||
{
|
||||
Eigen_matrix mat(nb_vert_, nb_vert_);
|
||||
|
||||
// cot values
|
||||
Eigen_matrix C = cot_entries();
|
||||
|
||||
for(face_descriptor f : faces(mesh_))
|
||||
{
|
||||
for(halfedge_descriptor hi : halfedges_around_face(halfedge(f, mesh_), mesh_))
|
||||
{
|
||||
vertex_descriptor v_source = source(hi, mesh_);
|
||||
vertex_descriptor v_target = target(hi, mesh_);
|
||||
auto i_source = vimap_[v_source];
|
||||
auto i_target = vimap_[v_target];
|
||||
|
||||
auto cot_val = C.coeff(i_source, i_target);
|
||||
// i,j
|
||||
mat.coeffRef(i_source, i_target) += cot_val; // why not (i,j)
|
||||
// j,i
|
||||
mat.coeffRef(i_target, i_source) += cot_val;
|
||||
//diagonal
|
||||
mat.coeffRef(i_source, i_source) -= cot_val;
|
||||
mat.coeffRef(i_target, i_target) -= cot_val;
|
||||
}
|
||||
}
|
||||
|
||||
return mat;
|
||||
}
|
||||
|
||||
|
||||
Eigen_matrix mass_matrix()
|
||||
{
|
||||
Eigen_matrix mat(nb_vert_, nb_vert_);
|
||||
for(face_descriptor f : faces(mesh_))
|
||||
{
|
||||
|
||||
double area = face_area(f, mesh_);
|
||||
|
||||
for(vertex_descriptor v : vertices_around_face(halfedge(f, mesh_), mesh_))
|
||||
{
|
||||
auto indx = vimap_[v];
|
||||
mat.coeffRef(indx, indx) += 2.0 * area;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
//export_eigen_matrix(mat, "data/D");
|
||||
|
||||
mat /= 12.0;
|
||||
|
||||
//export_eigen_matrix(mat, "data/D12");
|
||||
|
||||
return mat;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
void compute_coeff_matrix(Matrix& A, Eigen_matrix& L)
|
||||
{
|
||||
//Eigen_matrix L = get_stiffness_matrix();
|
||||
|
||||
//Eigen_matrix L = stiff_matrix();
|
||||
|
||||
export_eigen_matrix(L, "data/L");
|
||||
|
||||
//Eigen_matrix D = get_mass_matrix();
|
||||
Eigen_matrix D = mass_matrix();
|
||||
|
||||
export_eigen_matrix(D, "data/D");
|
||||
|
||||
|
||||
double delta = 0.001;
|
||||
|
||||
|
||||
//Eigen_matrix D(L.rows(), L.cols());
|
||||
//D.setIdentity();
|
||||
|
||||
|
||||
Eigen_matrix Ae = D - delta * L;
|
||||
|
||||
|
||||
fill_sparse_matrix(A, Ae);
|
||||
|
||||
A.assemble_matrix();
|
||||
|
||||
}
|
||||
|
||||
|
||||
void compute_rhs(Vector& Bx, Vector& By, Vector& Bz)
|
||||
{
|
||||
|
||||
Eigen_vector Xx(nb_vert_);
|
||||
Eigen_vector Xy(nb_vert_);
|
||||
Eigen_vector Xz(nb_vert_);
|
||||
|
||||
for(vertex_descriptor vi : vertices(mesh_))
|
||||
{
|
||||
int index = vimap_[vi];
|
||||
Point p = get(vpmap_, vi);
|
||||
Xx.coeffRef(index) = p.x();
|
||||
Xy.coeffRef(index) = p.y();
|
||||
Xz.coeffRef(index) = p.z();
|
||||
}
|
||||
|
||||
|
||||
//Eigen_matrix D = get_mass_matrix();
|
||||
Eigen_matrix D = mass_matrix();
|
||||
|
||||
|
||||
|
||||
//Eigen_matrix D(Bx.rows(), Bx.rows());
|
||||
//D.setIdentity();
|
||||
|
||||
|
||||
Eigen_vector Bxe = D * Xx;
|
||||
Eigen_vector Bye = D * Xy;
|
||||
Eigen_vector Bze = D * Xz;
|
||||
|
||||
|
||||
fill_sparse_vector(Bx, Bxe);
|
||||
fill_sparse_vector(By, Bye);
|
||||
fill_sparse_vector(Bz, Bze);
|
||||
|
||||
}
|
||||
|
||||
|
||||
// -------------- COPY TO SPARSE SYSTEM ---------- //
|
||||
void fill_sparse_matrix(Matrix&A, Eigen_matrix& Ae)
|
||||
{
|
||||
CGAL_assertion(A.row_dimension() == Ae.rows() &&
|
||||
A.column_dimension() == Ae.cols());
|
||||
|
||||
// this is not great
|
||||
for(std::size_t i=0; i<Ae.rows(); ++i)
|
||||
{
|
||||
for(std::size_t j=0; j<Ae.cols(); ++j)
|
||||
{
|
||||
A.set_coef(i, j, Ae.coeff(i, j), true);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void fill_sparse_vector(Vector& V, Eigen_vector& Ve)
|
||||
{
|
||||
CGAL_assertion(V.rows() == Ve.rows());
|
||||
for(std::size_t i= 0; i<Ve.rows(); ++i)
|
||||
{
|
||||
V[i] = Ve.coeff(i);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ---------------- CONSTRAINS ---------------- //
|
||||
void gather_constrained_vertices()
|
||||
{
|
||||
for(vertex_descriptor v : vertices(mesh_))
|
||||
{
|
||||
if(is_border(v, mesh_))
|
||||
{
|
||||
constrained_vertices_.insert(v);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void apply_constraints(Matrix& A)
|
||||
{ // to be called after gather_constrained_vertices()
|
||||
|
||||
typename std::unordered_set<vertex_descriptor>::iterator it;
|
||||
for(it = constrained_vertices_.begin(); it != constrained_vertices_.end(); ++it)
|
||||
{
|
||||
int i = get(vimap_, *it);
|
||||
A.set_coef(i, i, 1.0);
|
||||
|
||||
// also set all cols of the same row = 0 - not needed
|
||||
//for(std::size_t j = 0; j<A.column_dimension(); ++j)
|
||||
// A.set_coef(i, j, 0);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// -------------- EXTRACT : to move to debug file----------------- //
|
||||
void extract_matrix(Matrix& A)
|
||||
{
|
||||
std::ofstream out("data/mat.dat");
|
||||
for(int i=0; i < A.row_dimension(); ++i)
|
||||
{
|
||||
for(int j=0; j < A.column_dimension(); ++j)
|
||||
{
|
||||
NT val = A.get_coef(i, j);
|
||||
out<<val<<" ";
|
||||
}
|
||||
out<<nl;
|
||||
}
|
||||
out.close();
|
||||
}
|
||||
|
||||
void extract_vectors(Vector& Vx, Vector& Vy, Vector& Vz)
|
||||
{
|
||||
CGAL_assertion(Vx.rows() == Vy.rows());
|
||||
CGAL_assertion(Vx.rows() == Vz.rows());
|
||||
|
||||
std::ofstream out("data/vecs.dat");
|
||||
for(int j=0; j < Vx.dimension(); ++j)
|
||||
{
|
||||
out<<Vx[j]<<"\t"<<Vy[j]<<"\t"<<Vz[j]<<"\t"<<nl;
|
||||
}
|
||||
out.close();
|
||||
}
|
||||
|
||||
void export_eigen_matrix(Eigen_matrix& A, const char* filename)
|
||||
{
|
||||
|
||||
std::ofstream out(filename);
|
||||
for(auto i=0; i < A.rows(); ++i)
|
||||
{
|
||||
for(auto j=0; j < A.cols(); ++j)
|
||||
{
|
||||
NT val = A.coeff(i, j);
|
||||
out<<val<<" ";
|
||||
}
|
||||
out<<nl;
|
||||
}
|
||||
out.close();
|
||||
}
|
||||
|
||||
|
||||
// -------------- UPDATE MESH ----------------- //
|
||||
|
||||
Triangle triangle(face_descriptor f) const
|
||||
{
|
||||
halfedge_descriptor h = halfedge(f, mesh_);
|
||||
vertex_descriptor v1 = target(h, mesh_);
|
||||
vertex_descriptor v2 = target(next(h, mesh_), mesh_);
|
||||
vertex_descriptor v3 = target(next(next(h, mesh_), mesh_), mesh_);
|
||||
return Triangle(get(vpmap_, v1), get(vpmap_, v2), get(vpmap_, v3));
|
||||
}
|
||||
|
||||
void translate_centroid(Vector& Xx, Vector& Xy, Vector& Xz)
|
||||
{
|
||||
|
||||
std::vector<std::pair<Point, NT>> barycenters;
|
||||
|
||||
for(face_descriptor f : faces(mesh_))
|
||||
{
|
||||
Point tr_centroid = CGAL::centroid(triangle(f));
|
||||
barycenters.push_back(std::make_pair(tr_centroid, face_area(f, mesh_)));
|
||||
}
|
||||
|
||||
Point centroid = CGAL::barycenter(barycenters.begin(), barycenters.end());
|
||||
std::cout << "centroid= " << centroid << std::endl;
|
||||
|
||||
for(std::size_t i=0; i < Xx.rows(); ++i)
|
||||
{
|
||||
Xx[i] -= centroid.x();
|
||||
Xy[i] -= centroid.y();
|
||||
Xz[i] -= centroid.z();
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void normalize_area(Vector& Xx, Vector& Xy, Vector& Xz)
|
||||
{
|
||||
|
||||
NT surface_area = area(faces(mesh_), mesh_);
|
||||
|
||||
for(std::size_t i=0; i < Xx.rows(); ++i)
|
||||
{
|
||||
Xx[i] /= CGAL::sqrt(surface_area);
|
||||
Xy[i] /= CGAL::sqrt(surface_area);
|
||||
Xz[i] /= CGAL::sqrt(surface_area);
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
Xx /= surface_area;
|
||||
Xy /= surface_area;
|
||||
Xz /= surface_area;
|
||||
*/
|
||||
|
||||
}
|
||||
|
||||
|
||||
void update_map(Vector& Xx, Vector& Xy, Vector& Xz)
|
||||
{
|
||||
for (vertex_descriptor v : vertices(mesh_))
|
||||
{
|
||||
int index = get(vimap_, v);
|
||||
NT x_new = Xx[index];
|
||||
NT y_new = Xy[index];
|
||||
NT z_new = Xz[index];
|
||||
put(vpmap_, v, Point(x_new, y_new, z_new));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
} // PMP
|
||||
} // CGAL
|
||||
|
||||
|
||||
#endif // SMOOTHING_H
|
||||
|
|
@ -0,0 +1,70 @@
|
|||
#include <boost/graph/graph_traits.hpp>
|
||||
#include <boost/property_map/property_map.hpp>
|
||||
|
||||
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/implicit_shape_smoothing_impl.h>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
|
||||
namespace Polygon_mesh_processing {
|
||||
|
||||
|
||||
|
||||
template<typename PolygonMesh>
|
||||
void smooth_shape(PolygonMesh& mesh, int nb_iter)
|
||||
{
|
||||
|
||||
// VPmap type
|
||||
typedef typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type VertexPointMap;
|
||||
VertexPointMap vpmap = get(CGAL::vertex_point, mesh);
|
||||
|
||||
internal::Shape_smoother<PolygonMesh, VertexPointMap> smoother(mesh, vpmap);
|
||||
|
||||
Eigen::SparseMatrix<double> stiffness_matrix;
|
||||
stiffness_matrix = smoother.calc_stiff_matrix();
|
||||
|
||||
for(unsigned int t=0; t<nb_iter; ++t)
|
||||
{
|
||||
smoother.solve_system(stiffness_matrix);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
// demo API, undocumented
|
||||
template<typename PolygonMesh>
|
||||
void solve_mcf_system(PolygonMesh& mesh, int nb_iter, Eigen::SparseMatrix<double>& stiffness_matrix)
|
||||
{
|
||||
|
||||
// VPmap type
|
||||
typedef typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type VertexPointMap;
|
||||
VertexPointMap vpmap = get(CGAL::vertex_point, mesh);
|
||||
|
||||
internal::Shape_smoother<PolygonMesh, VertexPointMap> smoother(mesh, vpmap);
|
||||
|
||||
for(unsigned int t=0; t<nb_iter; ++t)
|
||||
{
|
||||
smoother.solve_system(stiffness_matrix);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
template<typename PolygonMesh>
|
||||
void setup_mcf_system(PolygonMesh& mesh, int nb_iter, Eigen::SparseMatrix<double>& stiffness_matrix)
|
||||
{
|
||||
|
||||
// VPmap type
|
||||
typedef typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type VertexPointMap;
|
||||
VertexPointMap vpmap = get(CGAL::vertex_point, mesh);
|
||||
|
||||
internal::Shape_smoother<PolygonMesh, VertexPointMap> smoother(mesh, vpmap);
|
||||
|
||||
stiffness_matrix = smoother.calc_stiff_matrix();
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
} //Polygon_mesh_processing
|
||||
|
||||
} //CGAL
|
||||
Loading…
Reference in New Issue