From 802102fa2e16f6caa9a57d9041089cef8d62ce0a Mon Sep 17 00:00:00 2001 From: Konstantinos Katrioplas Date: Tue, 13 Mar 2018 15:52:33 +0100 Subject: [PATCH] wip: testing perfomance --- .../internal/Smoothing/constraints_map.h | 66 +++++ .../internal/Smoothing/curvature_flow_impl.h | 93 +------ .../Smoothing/curvature_flow_new_impl.h | 263 +++++++----------- .../Polygon_mesh_processing/shape_smoothing.h | 110 +++++--- .../test_shape_smoothing.cpp | 37 ++- 5 files changed, 292 insertions(+), 277 deletions(-) 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 0005d5fbc50..24c80b25eba 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 @@ -33,6 +33,72 @@ namespace internal { } }; + + template> + 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::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::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 + 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 627f14986aa..87c8c619abd 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 @@ -38,70 +38,12 @@ #include +#include + namespace CGAL { namespace Polygon_mesh_processing { namespace internal { -/* -template> -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::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::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 -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; -}; -*/ template #include +#include namespace CGAL { namespace Polygon_mesh_processing { namespace internal { -template> -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::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::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 -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; -}; - template(), t.get<1>(), -time * t.get<2>()); + } + + + /* + std::vector diag_set(diagonal_.size(), false); + + for(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); + diag_set[t.get<0>()]=true; + } + } + + for(int idx = 0; idx < diagonal_.size(); ++idx) + { + if (!diag_set[idx]) + { + std::cout << "diag_NOT_set\n" ; + A.set_coef(idx, idx, diagonal_[idx], true); // A is empty + } + } + */ + + + A.assemble_matrix(); // does setFromTriplets and some compression + + } + + void calculate_D_diagonal(std::vector& diagonal) + { + diagonal.clear(); + //diagonal.resize(nb_vert_); + diagonal.assign(nb_vert_, 0); for(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; + } } } - D /= 12.0; + + for(int i = 0; i < diagonal.size(); ++i) + { + diagonal[i] /= 12.0; // constraints (= 1) as well ? + } + } - void compute_coefficient_matrix(Eigen_matrix& A, 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()); + 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::vector diagonal; - calculate_D_diagonal(diagonal); - - - - - - - //A = D - time * L; - } - - - void calculate_D_diagonal(std::vector& diagonal) - { - // to fill - } - - - void compute_rhs(Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz, - Eigen_matrix& D) - { for(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.set(index, diagonal_[index] * p.x()); + by.set(index, diagonal_[index] * p.y()); + bz.set(index, diagonal_[index] * p.z()); } - bx = D * bx; - by = D * by; - bz = D * bz; } // update mesh @@ -300,24 +257,14 @@ private: 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_)) { int index = get(vimap_, v); - NT x_new = Xx.coeffRef(index); - NT y_new = Xy.coeffRef(index); - NT z_new = Xz.coeffRef(index); + 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)); } } @@ -346,7 +293,11 @@ private: PolygonMesh& mesh_; VertexPointMap& vpmap_; std::size_t nb_vert_; - std::vector tripletList_; + + // 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_); 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 50ee6c8e37f..a4fc6b72ead 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 @@ -33,6 +33,9 @@ #include #include +#ifdef CGAL_PMP_SMOOTHING_VERBOSE +#include +#endif // new cgal traits #include @@ -221,7 +224,7 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con #ifdef CGAL_PMP_SMOOTHING_VERBOSE t.stop(); std::cout << " done ("<< t.time() <<" sec)." << std::endl; - std::cout << "Initializing..." << std::endl; + std::cout << "Initializing..."; t.reset(); t.start(); #endif @@ -231,31 +234,61 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con #ifdef CGAL_PMP_SMOOTHING_VERBOSE t.stop(); std::cout << " done ("<< t.time() <<" sec)." << std::endl; - std::cout << "#iter = " << nb_iterations << std::endl; - std::cout << "Solving linear system..." << 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 << "Shape smoothing done in "; + 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(const FaceRange& faces, PolygonMesh& pmesh, const double& time, +void smooth_along_curvature_flow_new(const FaceRange& faces, PolygonMesh& pmesh, const double& time, const NamedParameters& np) { using boost::choose_param; @@ -263,7 +296,7 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con #ifdef CGAL_PMP_SMOOTHING_VERBOSE CGAL::Timer t; - std::cout << "Smoothing parameters..."; + std::cout << "Initializing..."; std::cout.flush(); t.start(); #endif @@ -296,7 +329,9 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con // implicit scheme #if defined(CGAL_EIGEN3_ENABLED) #if EIGEN_VERSION_AT_LEAST(3,2,0) - typedef CGAL::Eigen_solver_traits<> Default_solver; + 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 @@ -320,63 +355,68 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con typedef typename Default_solver::Vector Eigen_vector; std::size_t n = static_cast(vertices(pmesh).size()); - Eigen_matrix A(n, n), stiffness_matrix(n, n), mass_matrix(n, n); + Eigen_matrix A(n, n); Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n); internal::Shape_smoother_new smoother(pmesh, vpmap, vcmap); + smoother.init_smoothing(faces); + #ifdef CGAL_PMP_SMOOTHING_VERBOSE t.stop(); std::cout << " done ("<< t.time() <<" sec)." << std::endl; - std::cout << "Initializing..." << std::endl; + std::cout << "calculate_stiffness_matrix_elements... "; t.reset(); t.start(); #endif - smoother.init_smoothing(faces); smoother.calculate_stiffness_matrix_elements(); #ifdef CGAL_PMP_SMOOTHING_VERBOSE t.stop(); std::cout << " done ("<< t.time() <<" sec)." << std::endl; - std::cout << "#iter = " << nb_iterations << std::endl; - std::cout << "Solving linear system..." << std::endl; - 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 + std::cout << "setup_system... "; + t.reset(); t.start(); + #endif + + smoother.setup_system(A, 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 ("<< t.time() <<" sec)." << std::endl; + #endif + } } - #ifdef CGAL_PMP_SMOOTHING_VERBOSE - t.stop(); - std::cout << "Shape smoothing 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()); -} - // demo helpers, undocumented template void setup_mcf_system(const FaceRange& faces, PolygonMesh& mesh, 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 af457336054..8d7466a7aa4 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 @@ -1,10 +1,11 @@ +#define CGAL_PMP_SMOOTHING_VERBOSE 1 +#define CGAL_PMP_SMOOTHING_DEBUG 1 #include #include #include #include #include - typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; typedef Kernel::Point_3 Point; typedef CGAL::Surface_mesh Mesh; @@ -140,6 +141,29 @@ 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) { std::ifstream input(filename); @@ -168,11 +192,12 @@ int main(int argc, char* argv[]) const char* filename_devil = "data/mannequin-devil.off"; const char* filename_pyramid = "data/simple_pyramid.off"; - test_demo_helpers(filename_devil); - test_curvature_flow_time_step(filename_devil); - test_curvature_flow(filename_pyramid); - test_implicit_constrained(filename_devil); - test_explicit_scheme(filename_devil); + //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_explicit_scheme(filename_devil); return 0; }