From e96d7207d4d368dba4d41d13da0939f5b78988a0 Mon Sep 17 00:00:00 2001 From: Konstantinos Katrioplas Date: Thu, 15 Mar 2018 13:14:58 +0100 Subject: [PATCH] testing all cases --- .../internal/Smoothing/constraints_map.h | 22 -- .../internal/Smoothing/curvature_flow_impl.h | 315 ++++++++---------- .../Smoothing/curvature_flow_new_impl.h | 305 ----------------- .../Polygon_mesh_processing/shape_smoothing.h | 211 +----------- .../test_shape_smoothing.cpp | 64 ++-- 5 files changed, 183 insertions(+), 734 deletions(-) delete mode 100644 Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_new_impl.h diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/constraints_map.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/constraints_map.h index 24c80b25eba..9bb6d93f2b4 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/constraints_map.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/constraints_map.h @@ -73,28 +73,6 @@ namespace internal { } }; - template - struct Incident_area - { - Incident_area(PolygonMesh& mesh) : pmesh(mesh){} - - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::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::null_face() ? 0 : face_area(f1, pmesh); - double A2 = f2 == boost::graph_traits::null_face() ? 0 : face_area(f2, pmesh); - return A1 + A2; - } - PolygonMesh& pmesh; - }; - - 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 8c386004f0f..4c38817fde3 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 @@ -18,28 +18,27 @@ // // Author(s) : Konstantinos Katrioplas (konst.katrioplas@gmail.com) -#ifndef CGAL_POLYGON_MESH_PROCESSING_CURVATURE_FLOW_IMPL_H -#define CGAL_POLYGON_MESH_PROCESSING_CURVATURE_FLOW_IMPL_H +#ifndef CURVATURE_FLOW_NEW_IMPL_H +#define CURVATURE_FLOW_NEW_IMPL_H + +#include +#include #include #include -#include -#include +#include #include #include - -#include -#include #include -#include - +#include +#if defined(CGAL_EIGEN3_ENABLED) +#include +#endif #include #include #include -#include - namespace CGAL { namespace Polygon_mesh_processing { namespace internal { @@ -47,6 +46,7 @@ namespace internal { template class Shape_smoother{ @@ -54,30 +54,26 @@ private: typedef typename GeomTraits::FT NT; typedef typename GeomTraits::Point_3 Point; - typedef typename GeomTraits::Triangle_3 Triangle; + typedef CGAL::Triple Triplet; typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits::face_descriptor face_descriptor; typedef typename boost::graph_traits::edge_descriptor edge_descriptor; - typedef typename boost::property_map::type IndexMap; - typedef typename Eigen::SparseMatrix Eigen_matrix; - typedef typename Eigen::VectorXd Eigen_vector; - - typedef typename Eigen::BiCGSTAB > Eigen_solver; - // typedef typename Eigen::SimplicialLDLT Eigen_solver; + // linear system + typedef typename SparseLinearSolver::Matrix Eigen_matrix; + typedef typename SparseLinearSolver::Vector Eigen_vector; public: Shape_smoother(PolygonMesh& mesh, - VertexPointMap& vpmap, - VertexConstraintMap& vcmap) : + VertexPointMap& vpmap, + VertexConstraintMap& vcmap) : mesh_(mesh), vpmap_(vpmap), vcmap_(vcmap), weight_calculator_(mesh, vpmap), - inc_areas_calculator_(mesh), nb_vert_(static_cast(vertices(mesh).size())) {} template @@ -86,216 +82,168 @@ public: check_face_range(face_range); } - void setup_system(Eigen_matrix& A, const Eigen_matrix& L, Eigen_matrix& D, + void setup_system(Eigen_matrix& A, Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz, + std::vector& stiffness_elements, const double& time) { - calc_mass_matrix(D); - - compute_coeff_matrix(A, L, D, time); - - compute_rhs(bx, by, bz, D); - + compute_coefficient_matrix(A, stiffness_elements, time); + compute_rhs(bx, by, bz); } void solve_system(const Eigen_matrix& A, Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz, const Eigen_vector& bx, const Eigen_vector& by, const Eigen_vector& bz) { + SparseLinearSolver solver; + NT D; - Eigen_solver solver; - solver.compute(A); - - - Xx = solver.solve(bx); - Xy = solver.solve(by); - Xz = solver.solve(bz); - - - - if(solver.info() != Eigen::Success) + // calls compute once to factorize with the preconditioner + if(!solver.factor(A, D)) { - std::cerr << "Not Solved!" << std::endl; + std::cerr << "Could not factorize linear system with preconditioner." << std::endl; return; } - - /* - std::cout << "Xx old:\n"; - for(int i = 0; i < Xx.rows(); ++i) + if(!solver.linear_solver(bx, Xx) || + !solver.linear_solver(by, Xy) || + !solver.linear_solver(bz, Xz) ) { - std::cout << Xx[i] << std::endl; + std::cerr << "Could not solve linear system." << std::endl; + return; } - std::cout <<"\n"; - */ } - void calc_stiff_matrix(Eigen_matrix& mat) + void calculate_stiffness_matrix_elements(std::vector& stiffness_elements) { - typedef Eigen::Triplet Triplet; - std::vector tripletList; - tripletList.reserve(8 * nb_vert_); - // todo: calculate exactly how many non zero entries there will be. + CGAL_assertion(stiffness_elements.empty()); + stiffness_elements.reserve(8 * nb_vert_); //estimation - for(face_descriptor f : frange_) + boost::unordered_map diag_coeff; + BOOST_FOREACH(face_descriptor f, frange_) { - for(halfedge_descriptor hi : halfedges_around_face(halfedge(f, mesh_), mesh_)) + BOOST_FOREACH(halfedge_descriptor hi, halfedges_around_face(halfedge(f, mesh_), mesh_)) { + halfedge_descriptor hi_opp = opposite(hi, mesh_); + if(!is_border(hi_opp, mesh_) && hi > hi_opp) continue; vertex_descriptor v_source = source(hi, mesh_); vertex_descriptor v_target = target(hi, mesh_); if(!is_constrained(v_source) && !is_constrained(v_target)) { - auto i_source = vimap_[v_source]; - auto i_target = vimap_[v_target]; + int i_source = vimap_[v_source]; + int i_target = vimap_[v_target]; NT Lij = weight_calculator_(hi); - tripletList.push_back(Triplet(i_source, i_target, Lij)); - tripletList.push_back(Triplet(i_target, i_source, Lij)); - tripletList.push_back(Triplet(i_source, i_source, -Lij)); - tripletList.push_back(Triplet(i_target, i_target, -Lij)); + if (!is_border(hi_opp, mesh_)) + Lij+= weight_calculator_(hi_opp); + stiffness_elements.push_back(Triplet(i_source, i_target, Lij)); + stiffness_elements.push_back(Triplet(i_target, i_source, Lij)); + diag_coeff.insert(std::make_pair(i_source,0)).first->second -= Lij; + diag_coeff.insert(std::make_pair(i_target,0)).first->second -= Lij; } } } - mat.setFromTriplets(tripletList.begin(), tripletList.end()); - - - + for(typename boost::unordered_map::iterator p = diag_coeff.begin(); + p != diag_coeff.end(); ++p) + { + stiffness_elements.push_back(Triplet(p->first, p->first, p->second)); + } } void update_mesh(Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz) { - update_map(Xx, Xy, Xz); - - // normalize area - //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); - //update_map(Xx, Xy, Xz); - //surface_area = area(faces(mesh_), mesh_); - //std::cout << "surface_area normalized= " << surface_area << nl; + 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)); + } } private: - void calc_mass_matrix(Eigen_matrix& D) + + // compute linear system + // ----------------------------------- + void compute_coefficient_matrix(Eigen_matrix& A, + std::vector& stiffness_elements, + const double& time) { - for(face_descriptor f : frange_) + CGAL_assertion(A.row_dimension() != 0); + CGAL_assertion(A.column_dimension() != 0); + CGAL_assertion(A.row_dimension() == nb_vert_); + CGAL_assertion(A.column_dimension() == nb_vert_); + + calculate_D_diagonal(diagonal_); + CGAL_assertion(diagonal_.size() == nb_vert_); + + // fill A = D - time * L + BOOST_FOREACH(Triplet t, stiffness_elements) + { + if (t.get<0>() != t.get<1>()) + A.set_coef(t.get<0>(), t.get<1>(), -time * t.get<2>(), true); + else + { + A.set_coef(t.get<0>(), t.get<1>(), diagonal_[t.get<0>()] -time * t.get<2>(), true); + } + } + A.assemble_matrix(); // does setFromTriplets and some compression + } + + void calculate_D_diagonal(std::vector& diagonal) + { + diagonal.clear(); + diagonal.assign(nb_vert_, 0); + std::vector constraints_flags(diagonal.size(), false); + BOOST_FOREACH(face_descriptor f, frange_) { double area = face_area(f, mesh_); for(vertex_descriptor v : vertices_around_face(halfedge(f, mesh_), mesh_)) { - auto idx = vimap_[v]; + int idx = vimap_[v]; if(!is_constrained(v)) - D.coeffRef(idx, idx) += 2.0 * area; + { + diagonal[idx] += 2.0 * area; + } else - D.coeffRef(idx, idx) = 1.0; + { + diagonal[idx] = 1.0; + constraints_flags[idx] = true; + } } } - D /= 12.0; - - /* - std::ofstream out_d("data/diag_eigen.cat"); - for(int i = 0; i < D.rows(); ++i) + for(int i = 0; i < diagonal.size(); ++i) { - out_d << D.coeffRef(i, i) << std::endl; + //if(!constraints_flags[i] == true) + diagonal[i] /= 12.0; } - */ - - - - } - void compute_coeff_matrix(Eigen_matrix& A, const Eigen_matrix& L, const Eigen_matrix& D, const double& time) + void compute_rhs(Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz) { - assert(A.rows() != 0); - assert(A.cols() != 0); - assert(A.rows() == L.rows()); - assert(A.cols() == L.cols()); - assert(A.rows() == D.rows()); - assert(A.cols() == D.cols()); - A = D - time * L; + CGAL_assertion(diagonal_.size() == nb_vert_); + CGAL_assertion(bx.size() == nb_vert_); + CGAL_assertion(by.size() == nb_vert_); + CGAL_assertion(bz.size() == nb_vert_); - /* - std::ofstream out("data/eigen_A.dat"); - for(auto j = 0 ; j < A.cols(); ++j) - { - for(auto i = 0; i < A.rows(); ++i) - { - - out << A.coeffRef(i, j) << " "; - } - out << std::endl; - } - */ - - - } - - void compute_rhs(Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz, - Eigen_matrix& D) - { - - for(vertex_descriptor vi : vrange_) + BOOST_FOREACH(vertex_descriptor vi, vrange_) { int index = vimap_[vi]; Point p = get(vpmap_, vi); - bx.coeffRef(index) = p.x(); - by.coeffRef(index) = p.y(); - bz.coeffRef(index) = p.z(); - } - bx = D * bx; - by = D * by; - bz = D * bz; - - - std::ofstream outbx("data/outb-eigen-x"); - std::ofstream outby("data/outb-eigen-y"); - std::ofstream outbz("data/outb-eigen-z"); - for(int i = 0 ; i < bx.rows(); ++i) - { - outbx << bx[i] << std::endl; - outby << by[i] << std::endl; - outbz << bz[i] << std::endl; - + bx.set(index, p.x()); + by.set(index, p.y()); + bz.set(index, p.z()); } - - - } - - // 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 normalize_area(Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz) - { - NT surface_area = area(faces(mesh_), mesh_); - Xx /= CGAL::sqrt(surface_area); - Xy /= CGAL::sqrt(surface_area); - Xz /= CGAL::sqrt(surface_area); - } - */ - - void update_map(Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz) - { - for (vertex_descriptor v : vertices(mesh_)) + // multiply D * b + for(int i = 0; i < nb_vert_; ++i) { - int index = get(vimap_, v); - NT x_new = Xx.coeffRef(index); - NT y_new = Xy.coeffRef(index); - NT z_new = Xz.coeffRef(index); - put(vpmap_, v, Point(x_new, y_new, z_new)); + bx[i] *= diagonal_[i]; + by[i] *= diagonal_[i]; + bz[i] *= diagonal_[i]; } } @@ -320,15 +268,19 @@ private: private: // data members // ------------ + std::size_t nb_vert_; PolygonMesh& mesh_; VertexPointMap& vpmap_; - std::size_t nb_vert_; + VertexConstraintMap& vcmap_; + IndexMap vimap_ = get(boost::vertex_index, mesh_); + + // linear system data + //std::vector stiffness_elements_; + std::vector diagonal_; // index of vector -> index of vimap_ + std::set frange_; std::set vrange_; - IndexMap vimap_ = get(boost::vertex_index, mesh_); - VertexConstraintMap vcmap_; Edge_cotangent_weight weight_calculator_; - Incident_area inc_areas_calculator_; }; @@ -338,4 +290,17 @@ private: } // CGAL -#endif // CGAL_POLYGON_MESH_PROCESSING_CURVATURE_FLOW_IMPL_H + + + + + + + + + + + + + +#endif // CURVATURE_FLOW_NEW_IMPL_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_new_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_new_impl.h deleted file mode 100644 index 1fba51faf98..00000000000 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_new_impl.h +++ /dev/null @@ -1,305 +0,0 @@ -// Copyright (c) 2017 GeometryFactory (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// You can redistribute it and/or modify it under the terms of the GNU -// General Public License as published by the Free Software Foundation, -// either version 3 of the License, or (at your option) any later version. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL$ -// $Id$ -// -// -// Author(s) : Konstantinos Katrioplas (konst.katrioplas@gmail.com) - - -#ifndef CURVATURE_FLOW_NEW_IMPL_H -#define CURVATURE_FLOW_NEW_IMPL_H - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#if defined(CGAL_EIGEN3_ENABLED) -#include -#endif -#include -#include -#include - - -namespace CGAL { -namespace Polygon_mesh_processing { -namespace internal { - -template -class Shape_smoother_new{ - -private: - - typedef typename GeomTraits::FT NT; - typedef typename GeomTraits::Point_3 Point; - typedef CGAL::Triple Triplet; - - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; - typedef typename boost::graph_traits::edge_descriptor edge_descriptor; - typedef typename boost::property_map::type IndexMap; - - // linear system - typedef typename SparseLinearSolver::Matrix Eigen_matrix; - typedef typename SparseLinearSolver::Vector Eigen_vector; - -public: - Shape_smoother_new(PolygonMesh& mesh, - VertexPointMap& vpmap, - VertexConstraintMap& vcmap) : - mesh_(mesh), - vpmap_(vpmap), - vcmap_(vcmap), - weight_calculator_(mesh, vpmap), - nb_vert_(static_cast(vertices(mesh).size())) {} - - template - void init_smoothing(const FaceRange& face_range) - { - check_face_range(face_range); - } - - void setup_system(Eigen_matrix& A, - Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz, - std::vector& stiffness_elements, - const double& time) - { - compute_coefficient_matrix(A, stiffness_elements, time); - compute_rhs(bx, by, bz); - } - - void solve_system(const Eigen_matrix& A, - Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz, - const Eigen_vector& bx, const Eigen_vector& by, const Eigen_vector& bz) - { - SparseLinearSolver solver; - NT D; - - // calls compute once to factorize with the preconditioner - if(!solver.factor(A, D)) - { - std::cerr << "Could not factorize lienar system with preconditioner." << std::endl; - } - - if(!solver.linear_solver(bx, Xx) || - !solver.linear_solver(by, Xy) || - !solver.linear_solver(bz, Xz) ) - { - std::cerr << "Could not solve linear system." << std::endl; - return; - } - } - - void calculate_stiffness_matrix_elements(std::vector& stiffness_elements) - { - CGAL_assertion(stiffness_elements.empty()); - stiffness_elements.reserve(8 * nb_vert_); //estimation - - boost::unordered_map diag_coeff; - BOOST_FOREACH(face_descriptor f, frange_) - { - BOOST_FOREACH(halfedge_descriptor hi, halfedges_around_face(halfedge(f, mesh_), mesh_)) - { - halfedge_descriptor hi_opp = opposite(hi, mesh_); - if(!is_border(hi_opp, mesh_) && hi > hi_opp) continue; - vertex_descriptor v_source = source(hi, mesh_); - vertex_descriptor v_target = target(hi, mesh_); - - if(!is_constrained(v_source) && !is_constrained(v_target)) - { - int i_source = vimap_[v_source]; - int i_target = vimap_[v_target]; - NT Lij = weight_calculator_(hi); - if (!is_border(hi_opp, mesh_)) - Lij+= weight_calculator_(hi_opp); - stiffness_elements.push_back(Triplet(i_source, i_target, Lij)); - stiffness_elements.push_back(Triplet(i_target, i_source, Lij)); - diag_coeff.insert(std::make_pair(i_source,0)).first->second -= Lij; - diag_coeff.insert(std::make_pair(i_target,0)).first->second -= Lij; - } - } - } - - for(typename boost::unordered_map::iterator p = diag_coeff.begin(); - p != diag_coeff.end(); ++p) - { - stiffness_elements.push_back(Triplet(p->first, p->first, p->second)); - } - } - - void update_mesh(Eigen_vector& Xx, Eigen_vector& Xy, Eigen_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)); - } - } - -private: - - // compute linear system - // ----------------------------------- - void compute_coefficient_matrix(Eigen_matrix& A, - std::vector& stiffness_elements, - const double& time) - { - CGAL_assertion(A.row_dimension() != 0); - CGAL_assertion(A.column_dimension() != 0); - CGAL_assertion(A.row_dimension() == nb_vert_); - CGAL_assertion(A.column_dimension() == nb_vert_); - - calculate_D_diagonal(diagonal_); - CGAL_assertion(diagonal_.size() == nb_vert_); - - // fill A = D - time * L - BOOST_FOREACH(Triplet t, stiffness_elements) - { - if (t.get<0>() != t.get<1>()) - A.set_coef(t.get<0>(), t.get<1>(), -time * t.get<2>(), true); - else - { - A.set_coef(t.get<0>(), t.get<1>(), diagonal_[t.get<0>()] -time * t.get<2>(), true); - } - } - A.assemble_matrix(); // does setFromTriplets and some compression - } - - void calculate_D_diagonal(std::vector& diagonal) - { - diagonal.clear(); - diagonal.assign(nb_vert_, 0); - std::vector constraints_flags(diagonal.size(), false); - BOOST_FOREACH(face_descriptor f, frange_) - { - double area = face_area(f, mesh_); - for(vertex_descriptor v : vertices_around_face(halfedge(f, mesh_), mesh_)) - { - int idx = vimap_[v]; - if(!is_constrained(v)) - { - diagonal[idx] += 2.0 * area; - } - else - { - diagonal[idx] = 1.0; - constraints_flags[idx] = true; - } - } - } - - for(int i = 0; i < diagonal.size(); ++i) - { - if(!constraints_flags[i] == true) - diagonal[i] /= 12.0; - } - } - - void compute_rhs(Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz) - { - CGAL_assertion(diagonal_.size() == nb_vert_); - CGAL_assertion(bx.size() == nb_vert_); - CGAL_assertion(by.size() == nb_vert_); - CGAL_assertion(bz.size() == nb_vert_); - - BOOST_FOREACH(vertex_descriptor vi, vrange_) - { - int index = vimap_[vi]; - Point p = get(vpmap_, vi); - bx.set(index, p.x()); - by.set(index, p.y()); - bz.set(index, p.z()); - } - - // multiply D * b - for(int i = 0; i < nb_vert_; ++i) - { - bx[i] *= diagonal_[i]; - by[i] *= diagonal_[i]; - bz[i] *= diagonal_[i]; - } - } - - // handling constrains - // ----------------------------------- - bool is_constrained(const vertex_descriptor& v) - { - return get(vcmap_, v); - } - - template - void check_face_range(const FaceRange& face_range) - { - BOOST_FOREACH(face_descriptor f, face_range) - { - frange_.insert(f); - BOOST_FOREACH(vertex_descriptor v, vertices_around_face(halfedge(f, mesh_), mesh_)) - vrange_.insert(v); - } - } - - private: - // data members - // ------------ - std::size_t nb_vert_; - PolygonMesh& mesh_; - VertexPointMap& vpmap_; - VertexConstraintMap& vcmap_; - IndexMap vimap_ = get(boost::vertex_index, mesh_); - - // linear system data - //std::vector stiffness_elements_; - std::vector diagonal_; // index of vector -> index of vimap_ - - std::set frange_; - std::set vrange_; - Edge_cotangent_weight weight_calculator_; - -}; - - -} // internal -} // PMP -} // CGAL - - - - - - - - - - - - - - - -#endif // CURVATURE_FLOW_NEW_IMPL_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_smoothing.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_smoothing.h index c336ed18ada..1b0b5794c31 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_smoothing.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_smoothing.h @@ -25,7 +25,6 @@ #include #if defined(CGAL_EIGEN3_ENABLED) -#include #include #endif @@ -37,8 +36,6 @@ #ifdef CGAL_PMP_SMOOTHING_VERBOSE #include #endif -// new cgal traits -#include namespace CGAL { namespace Polygon_mesh_processing { @@ -179,122 +176,6 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con using boost::choose_param; using boost::get_param; - #ifdef CGAL_PMP_SMOOTHING_VERBOSE - CGAL::Timer t; - std::cout << "Smoothing parameters..."; - std::cout.flush(); - t.start(); - #endif - - 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, pmesh)); - - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::lookup_named_param_def < - internal_np::vertex_is_constrained_t, - NamedParameters, - internal::Constrained_vertices_map - > ::type VCMap; - VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained), - internal::Constrained_vertices_map()); - - 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); - - if(use_explicit_scheme) - { - smooth_curvature_flow_explicit(faces, pmesh, parameters::number_of_iterations(nb_iterations)); - } - else - { - // implicit scheme - typedef typename Eigen::VectorXd Eigen_vector; - typedef typename Eigen::SparseMatrix Eigen_matrix; - - std::size_t n = static_cast(vertices(pmesh).size()); - Eigen_matrix A(n, n), stiffness_matrix(n, n), mass_matrix(n, n); - Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n); - - internal::Shape_smoother - smoother(pmesh, vpmap, vcmap); - - #ifdef CGAL_PMP_SMOOTHING_VERBOSE - t.stop(); - std::cout << " done ("<< t.time() <<" sec)." << std::endl; - std::cout << "Initializing..."; - t.reset(); t.start(); - #endif - - smoother.init_smoothing(faces); - smoother.calc_stiff_matrix(stiffness_matrix); - - #ifdef CGAL_PMP_SMOOTHING_VERBOSE - t.stop(); - std::cout << " done ("<< t.time() <<" sec)." << std::endl; - std::cout << "Setup system..."; - t.reset(); t.start(); - #endif - - for(std::size_t iter = 0; iter < nb_iterations; ++iter) - { - smoother.setup_system(A, stiffness_matrix, mass_matrix, bx, by, bz, time); - - #ifdef CGAL_PMP_SMOOTHING_VERBOSE - t.stop(); - std::cout << " done ("<< t.time() <<" sec)." << std::endl; - std::cout << "Solve system..."; - t.reset(); t.start(); - #endif - - smoother.solve_system(A, Xx, Xy, Xz, bx, by, bz); - - #ifdef CGAL_PMP_SMOOTHING_VERBOSE - t.stop(); - std::cout << " done ("<< t.time() <<" sec)." << std::endl; - std::cout << "update_mesh..."; - t.reset(); t.start(); - #endif - - - smoother.update_mesh(Xx, Xy, Xz); - } - } - - #ifdef CGAL_PMP_SMOOTHING_VERBOSE - t.stop(); - std::cout << " done in "; - std::cout << t.time() << " sec." << std::endl; - std::cout< -void smooth_along_curvature_flow(PolygonMesh& pmesh, const double& time, - const NamedParameters& np) -{ - smooth_along_curvature_flow(faces(pmesh), pmesh, time, np); -} - -template -void smooth_along_curvature_flow(PolygonMesh& pmesh, const double& time) -{ - smooth_along_curvature_flow(faces(pmesh), pmesh, time, - parameters::all_default()); -} - - -// new with cgal solver traits -template -void smooth_along_curvature_flow_new(const FaceRange& faces, PolygonMesh& pmesh, const double& time, - const NamedParameters& np) -{ - using boost::choose_param; - using boost::get_param; - #ifdef CGAL_PMP_SMOOTHING_VERBOSE CGAL::Timer t; std::cout << "Initializing..."; @@ -360,7 +241,7 @@ void smooth_along_curvature_flow_new(const FaceRange& faces, PolygonMesh& pmesh, Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n); std::vector > stiffness; - internal::Shape_smoother_new + internal::Shape_smoother smoother(pmesh, vpmap, vcmap); smoother.init_smoothing(faces); @@ -416,94 +297,20 @@ void smooth_along_curvature_flow_new(const FaceRange& faces, PolygonMesh& pmesh, } template -void smooth_along_curvature_flow_new(PolygonMesh& pmesh, const double& time, - const NamedParameters& np) +void smooth_along_curvature_flow(PolygonMesh& pmesh, const double& time, + const NamedParameters& np) { - smooth_along_curvature_flow_new(faces(pmesh), pmesh, time, np); + smooth_along_curvature_flow(faces(pmesh), pmesh, time, np); } template -void smooth_along_curvature_flow_new(PolygonMesh& pmesh, const double& time) +void smooth_along_curvature_flow(PolygonMesh& pmesh, const double& time) { - smooth_along_curvature_flow_new(faces(pmesh), pmesh, time, - parameters::all_default()); + smooth_along_curvature_flow(faces(pmesh), pmesh, time, + parameters::all_default()); } // API for Polyhedron demo plugin -template -void setup_mcf_system(const FaceRange& faces, PolygonMesh& mesh, - Eigen::SparseMatrix& stiffness_matrix, const NamedParameters& np) -{ - using boost::choose_param; - using boost::get_param; - - 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, mesh)); - - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::lookup_named_param_def < - internal_np::vertex_is_constrained_t, - NamedParameters, - internal::Constrained_vertices_map - > ::type VCMap; - VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained), - internal::Constrained_vertices_map()); - - std::size_t n = static_cast(vertices(mesh).size()); - stiffness_matrix.resize(n, n); - - internal::Shape_smoother - smoother(mesh, vpmap, vcmap); - smoother.init_smoothing(faces); - smoother.calc_stiff_matrix(stiffness_matrix); -} - -template -void solve_mcf_system(const FaceRange& faces, PolygonMesh& mesh, const double& time, - Eigen::SparseMatrix& stiffness_matrix, const NamedParameters& np) -{ - - using boost::choose_param; - using boost::get_param; - - //typedef typename GetGeomTraits::type GeomTraits; - //typedef typename boost::property_map::type VertexPointMap; - //VertexPointMap vpmap = get(CGAL::vertex_point, mesh); - - 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, mesh)); - - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::lookup_named_param_def < - internal_np::vertex_is_constrained_t, - NamedParameters, - internal::Constrained_vertices_map - > ::type VCMap; - VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained), - internal::Constrained_vertices_map()); - - typedef typename Eigen::VectorXd Eigen_vector; - typedef typename Eigen::SparseMatrix Eigen_matrix; - - std::size_t n = static_cast(vertices(mesh).size()); - - Eigen_matrix A(n, n), mass_matrix(n, n); - Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n); - - internal::Shape_smoother - smoother(mesh, vpmap, vcmap); - smoother.init_smoothing(faces); - smoother.setup_system(A, stiffness_matrix, mass_matrix, bx, by, bz, time); - smoother.solve_system(A, Xx, Xy, Xz, bx, by, bz); - smoother.update_mesh(Xx, Xy, Xz); -} - template void solve_mcf(const FaceRange& faces, PolygonMesh& mesh, const double& time, std::vector >& stiffness, bool compute_stiffness, @@ -565,14 +372,14 @@ void solve_mcf(const FaceRange& faces, PolygonMesh& mesh, const double& time, if(compute_stiffness) { - internal::Shape_smoother_new + internal::Shape_smoother smoother(mesh, vpmap, vcmap); smoother.init_smoothing(faces); smoother.calculate_stiffness_matrix_elements(stiffness); } else { - internal::Shape_smoother_new + internal::Shape_smoother smoother(mesh, vpmap, vcmap); smoother.init_smoothing(faces); smoother.setup_system(A, bx, by, bz, stiffness, time); 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 8d7466a7aa4..bb517139ab9 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 @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -47,6 +48,10 @@ public: void test_implicit_constrained(const char* filename) { + #ifdef CGAL_PMP_SMOOTHING_VERBOSE + std::cout << "-- test_implicit_constrained --" << std::endl; + #endif + std::ifstream input(filename); Mesh mesh; input >> mesh; @@ -70,7 +75,7 @@ void test_implicit_constrained(const char* filename) const double time_step = 1.0; CGAL::Polygon_mesh_processing::smooth_along_curvature_flow(mesh, time_step, CGAL::Polygon_mesh_processing::parameters::vertex_is_constrained_map(vcmap). - number_of_iterations(5)); + number_of_iterations(1)); #ifdef CGAL_PMP_SMOOTHING_VERBOSE std::ofstream out("data/output_implicit_constrained.off"); out << mesh; @@ -80,6 +85,10 @@ void test_implicit_constrained(const char* filename) void test_explicit_scheme(const char* filename) { + #ifdef CGAL_PMP_SMOOTHING_VERBOSE + std::cout << "-- test_explicit_scheme --" << std::endl; + #endif + std::ifstream input(filename); Mesh mesh; input >> mesh; @@ -103,6 +112,10 @@ void test_explicit_scheme(const char* filename) void test_curvature_flow_time_step(const char* filename) { + #ifdef CGAL_PMP_SMOOTHING_VERBOSE + std::cout << "-- test_curvature_flow_time_step --" << std::endl; + #endif + std::ifstream input(filename); Mesh mesh; input >> mesh; @@ -123,6 +136,10 @@ void test_curvature_flow_time_step(const char* filename) void test_curvature_flow(const char* filename) { + #ifdef CGAL_PMP_SMOOTHING_VERBOSE + std::cout << "-- test_curvature_flow --" << std::endl; + #endif + std::ifstream input(filename); Mesh mesh; input >> mesh; @@ -141,31 +158,12 @@ void test_curvature_flow(const char* filename) #endif } -void test_curvature_flow_new(const char* filename) -{ - std::ifstream input(filename); - Mesh mesh; - input >> mesh; - input.close(); - - boost::property_map::type vpmap = - get(CGAL::vertex_point, mesh); - - const double time_step = 1.0; - CGAL::Polygon_mesh_processing::smooth_along_curvature_flow_new(faces(mesh), mesh, time_step, - CGAL::Polygon_mesh_processing::parameters::all_default()); - - #ifdef CGAL_PMP_SMOOTHING_VERBOSE - std::ofstream out("data/output_precision_pyramid_new.off"); - out << mesh; - out.close(); - #endif -} - - - void test_demo_helpers(const char* filename) { + #ifdef CGAL_PMP_SMOOTHING_VERBOSE + std::cout << "-- test_demo_helpers --" << std::endl; + #endif + std::ifstream input(filename); Mesh mesh; input >> mesh; @@ -175,9 +173,16 @@ void test_demo_helpers(const char* filename) get(CGAL::vertex_point, mesh); const double time_step = 1e-2; - Eigen::SparseMatrix stiff_matrix; - CGAL::Polygon_mesh_processing::setup_mcf_system(faces(mesh), mesh, stiff_matrix, CGAL::Polygon_mesh_processing::parameters::all_default()); - CGAL::Polygon_mesh_processing::solve_mcf_system(faces(mesh), mesh, time_step, stiff_matrix, CGAL::Polygon_mesh_processing::parameters::all_default()); + std::vector > stiffness; + + bool compute_stiffness = true; + CGAL::Polygon_mesh_processing::solve_mcf(faces(mesh), mesh, time_step, + stiffness, compute_stiffness, + CGAL::Polygon_mesh_processing::parameters::all_default()); + compute_stiffness = false; + CGAL::Polygon_mesh_processing::solve_mcf(faces(mesh), mesh, time_step, + stiffness, compute_stiffness, + CGAL::Polygon_mesh_processing::parameters::all_default()); #ifdef CGAL_PMP_SMOOTHING_VERBOSE std::ofstream out("data/output_devil_demo_helpers.off"); @@ -194,9 +199,8 @@ int main(int argc, char* argv[]) //test_demo_helpers(filename_devil); //test_curvature_flow_time_step(filename_devil); - test_curvature_flow(filename_devil); - test_curvature_flow_new(filename_devil); - //test_implicit_constrained(filename_devil); + //test_curvature_flow(filename_pyramid); + test_implicit_constrained(filename_devil); //test_explicit_scheme(filename_devil); return 0;