From af7fa2ac51f906fc726452b5a8b51baf1286f019 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 10 Mar 2015 10:44:57 +0100 Subject: [PATCH 1/6] Add LSCM as parametrization method --- Polyhedron/demo/Polyhedron/MainWindow.ui | 12 +++++++++--- ...olyhedron_demo_parameterization_plugin.cpp | 19 +++++++++++++++++-- 2 files changed, 26 insertions(+), 5 deletions(-) diff --git a/Polyhedron/demo/Polyhedron/MainWindow.ui b/Polyhedron/demo/Polyhedron/MainWindow.ui index 1e183f30a21..581632e84ed 100644 --- a/Polyhedron/demo/Polyhedron/MainWindow.ui +++ b/Polyhedron/demo/Polyhedron/MainWindow.ui @@ -37,7 +37,7 @@ 0 0 978 - 21 + 26 @@ -94,6 +94,7 @@ + @@ -294,8 +295,8 @@ 0 0 - 545 - 178 + 534 + 174 @@ -646,6 +647,11 @@ &Preferences + + + Least square conformal maps + + diff --git a/Polyhedron/demo/Polyhedron/Polyhedron_demo_parameterization_plugin.cpp b/Polyhedron/demo/Polyhedron/Polyhedron_demo_parameterization_plugin.cpp index 8b3d53e2718..880c0e57019 100644 --- a/Polyhedron/demo/Polyhedron/Polyhedron_demo_parameterization_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Polyhedron_demo_parameterization_plugin.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include @@ -34,7 +35,8 @@ public: // used by Polyhedron_demo_plugin_helper QStringList actionsNames() const { return QStringList() << "actionMVC" - << "actionDCP"; + << "actionDCP" + << "actionLSC"; } bool applicable(QAction*) const { @@ -44,9 +46,10 @@ public: public slots: void on_actionMVC_triggered(); void on_actionDCP_triggered(); + void on_actionLSC_triggered(); protected: - enum Parameterization_method { PARAM_MVC, PARAM_DCP }; + enum Parameterization_method { PARAM_MVC, PARAM_DCP, PARAM_LSC }; void parameterize(Parameterization_method method); }; // end Polyhedron_demo_parameterization_plugin @@ -90,6 +93,13 @@ void Polyhedron_demo_parameterization_plugin::parameterize(const Parameterizatio typedef CGAL::Discrete_conformal_map_parameterizer_3 Parameterizer; Parameterizer::Error_code err = CGAL::parameterize(adaptor,Parameterizer()); success = err == Parameterizer::OK; + } + case PARAM_LSC: + { + std::cout << "Parameterize (LSC)..."; + typedef CGAL::LSCM_parameterizer_3 Parameterizer; + Parameterizer::Error_code err = CGAL::parameterize(adaptor,Parameterizer()); + success = err == Parameterizer::OK; } } @@ -145,6 +155,11 @@ void Polyhedron_demo_parameterization_plugin::on_actionDCP_triggered() parameterize(PARAM_DCP); } +void Polyhedron_demo_parameterization_plugin::on_actionLSC_triggered() +{ + std::cerr << "LSC..."; + parameterize(PARAM_LSC); +} Q_EXPORT_PLUGIN2(Polyhedron_demo_parameterization_plugin, Polyhedron_demo_parameterization_plugin) #include "Polyhedron_demo_parameterization_plugin.moc" From 18d3dda0a6efa00f76101b84876a283403b799d8 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 10 Mar 2015 10:48:22 +0100 Subject: [PATCH 2/6] Use a faster solver in the example --- .../Eigen_parameterization.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Eigen_parameterization.cpp b/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Eigen_parameterization.cpp index 7c97b3f8823..0d0300fbee9 100644 --- a/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Eigen_parameterization.cpp +++ b/Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Eigen_parameterization.cpp @@ -3,7 +3,7 @@ #include #include #include - +#include #include #include @@ -17,7 +17,7 @@ typedef CGAL::Simple_cartesian Kernel; typedef CGAL::Polyhedron_3 Polyhedron; - +typedef CGAL::Timer Timer; // ---------------------------------------------------------------------------- // main() @@ -65,7 +65,8 @@ int main(int argc, char * argv[]) typedef CGAL::Parameterization_polyhedron_adaptor_3 Parameterization_polyhedron_adaptor; - + Timer t; + t.start(); Parameterization_polyhedron_adaptor mesh_adaptor(mesh); //*************************************** @@ -77,16 +78,17 @@ int main(int argc, char * argv[]) typedef CGAL::Circular_border_arc_length_parameterizer_3 Border_parameterizer; // Eigen solver - typedef CGAL::Eigen_solver_traits<> Solver; + typedef CGAL::Eigen_solver_traits::EigenType, Eigen::IncompleteLUT< double > > > Solver; // Floater Mean Value Coordinates parameterization // (circular border) with Eigen solver typedef CGAL::Mean_value_coordinates_parameterizer_3 + Solver> Parameterizer; Parameterizer::Error_code err = CGAL::parameterize(mesh_adaptor, Parameterizer()); + t.stop(); switch(err) { case Parameterizer::OK: // Success @@ -119,6 +121,6 @@ int main(int argc, char * argv[]) double v = mesh_adaptor.info(pVertex->halfedge())->uv().y(); std::cout << "(u,v) = (" << u << "," << v << ")" << std::endl; } - + std::cerr << t.time() << "sec." << std::endl; return EXIT_SUCCESS; } From db4a974f640fb06f42e52c4ef0f5536ed9aecbd5 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Fri, 7 Aug 2015 11:55:00 +0200 Subject: [PATCH 3/6] Eigen interface for LSCM --- .../include/CGAL/Eigen_solver_traits.h | 5 ++ .../include/CGAL/LSCM_parameterizer_3.h | 88 +++++++++++-------- 2 files changed, 54 insertions(+), 39 deletions(-) diff --git a/Solver_interface/include/CGAL/Eigen_solver_traits.h b/Solver_interface/include/CGAL/Eigen_solver_traits.h index e641c121b89..9e065c65558 100644 --- a/Solver_interface/include/CGAL/Eigen_solver_traits.h +++ b/Solver_interface/include/CGAL/Eigen_solver_traits.h @@ -53,6 +53,11 @@ namespace internal { typedef Eigen_sparse_symmetric_matrix type; }; + template + struct Get_eigen_matrix< ::Eigen::LeastSquaresConjugateGradient,FT>{ + typedef Eigen_sparse_matrix type; + }; + template struct Get_eigen_matrix< ::Eigen::SimplicialCholesky,FT>{ typedef Eigen_sparse_symmetric_matrix type; diff --git a/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h index ef820d9b813..47d42c98f34 100644 --- a/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h @@ -24,14 +24,22 @@ #include #include + +#ifdef CGAL_EIGEN3_ENABLED +#include +#endif + #include + #include #include #include #include #include +#define DEBUG_TRACE + /// \file LSCM_parameterizer_3.h namespace CGAL { @@ -136,8 +144,12 @@ private: typedef typename Sparse_LA::Vector Vector; typedef typename Sparse_LA::Matrix Matrix; - typedef typename OpenNL::LinearSolver - LeastSquaresSolver ; + typedef Eigen_solver_traits::EigenType> > + LeastSquaresSolver ; + + typedef typename LeastSquaresSolver::Matrix LSS_Matrix; + typedef typename LeastSquaresSolver::Vector LSS_Vector; + typedef typename LeastSquaresSolver::NT LSS_Scalar; // Public operations public: @@ -173,7 +185,8 @@ private: /// \pre Vertices must be indexed. /// \pre X and B must be allocated and empty. /// \pre At least 2 border vertices must be parameterized. - void initialize_system_from_mesh_border(LeastSquaresSolver& solver, + void initialize_system_from_mesh_border(LSS_Matrix& A, + LSS_Vector& B, const Adaptor& mesh) ; /// Utility for setup_triangle_relations(): @@ -185,7 +198,7 @@ private: /// Create two lines in the linear system per triangle (one for u, one for v). /// /// \pre vertices must be indexed. - Error_code setup_triangle_relations(LeastSquaresSolver& solver, + Error_code setup_triangle_relations(LSS_Matrix& A, LSS_Vector& B, const Adaptor& mesh, Facet_const_handle facet) ; @@ -288,34 +301,39 @@ parameterize(Adaptor& mesh) // Create sparse linear system "A*X = B" of size 2*nbVertices x 2*nbVertices // (in fact, we need only 2 lines per triangle x 1 column per vertex) - LeastSquaresSolver solver(2*nbVertices); - solver.set_least_squares(true) ; - + LSS_Matrix A (2 * nbVertices, 2 * nbVertices); + LSS_Vector B (2 * nbVertices); + // Initialize the "A*X = B" linear system after // (at least two) border vertices parameterization - initialize_system_from_mesh_border(solver, mesh); + initialize_system_from_mesh_border(A, B, mesh); // Fill the matrix for the other vertices - solver.begin_system() ; for (Facet_iterator facetIt = mesh.mesh_facets_begin(); facetIt != mesh.mesh_facets_end(); facetIt++) { // Create two lines in the linear system per triangle (one for u, one for v) - status = setup_triangle_relations(solver, mesh, facetIt); + status = setup_triangle_relations(A, B, mesh, facetIt); if (status != Base::OK) return status; } - solver.end_system() ; + + A.assemble_matrix(); + #ifdef DEBUG_TRACE std::cerr << " matrix filling (" << 2*mesh.count_mesh_facets() << " x " << nbVertices << "): " << timer.time() << " seconds." << std::endl; timer.reset(); #endif - + // Solve the "A*X = B" linear system in the least squares sense - if ( ! solver.solve() ) - status = Base::ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; + LSS_Vector X (2 * nbVertices); + LSS_Scalar D; + LeastSquaresSolver solver; + if (!(solver.linear_solver (A, B, X, D))) + status = Base::ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; + #ifdef DEBUG_TRACE std::cerr << " solving linear system: " << timer.time() << " seconds." << std::endl; @@ -323,7 +341,7 @@ parameterize(Adaptor& mesh) #endif if (status != Base::OK) return status; - + return status; // Copy X coordinates into the (u,v) pair of each vertex set_mesh_uv_from_system(mesh, solver); #ifdef DEBUG_TRACE @@ -396,7 +414,7 @@ check_parameterize_preconditions(Adaptor& mesh) template inline void LSCM_parameterizer_3:: -initialize_system_from_mesh_border(LeastSquaresSolver& solver, +initialize_system_from_mesh_border(LSS_Matrix& A, LSS_Vector& B, const Adaptor& mesh) { for (Vertex_const_iterator it = mesh.mesh_vertices_begin(); @@ -412,14 +430,13 @@ initialize_system_from_mesh_border(LeastSquaresSolver& solver, // Write (u,v) in X (meaningless if vertex is not parameterized) // Note : 2*index --> u // 2*index + 1 --> v - solver.variable(2*index ).set_value(uv.x()) ; - solver.variable(2*index + 1).set_value(uv.y()) ; // Copy (u,v) in B if vertex is parameterized - if (mesh.is_vertex_parameterized(it)) { - solver.variable(2*index ).lock() ; - solver.variable(2*index + 1).lock() ; - } + if (mesh.is_vertex_parameterized(it)) + { + B.set (2 * index, uv.x ()); + B.set (2 * index + 1, uv.y ()); + } } } @@ -474,7 +491,7 @@ template inline typename LSCM_parameterizer_3::Error_code LSCM_parameterizer_3:: -setup_triangle_relations(LeastSquaresSolver& solver, +setup_triangle_relations(LSS_Matrix& A, LSS_Vector& B, const Adaptor& mesh, Facet_const_handle facet) { @@ -537,24 +554,15 @@ setup_triangle_relations(LeastSquaresSolver& solver, // // Real part // Note: b = 0 - solver.begin_row() ; - solver.add_coefficient(u0_id, -a+c) ; - solver.add_coefficient(v0_id, b-d) ; - solver.add_coefficient(u1_id, -c) ; - solver.add_coefficient(v1_id, d) ; - solver.add_coefficient(u2_id, a) ; - solver.end_row() ; + + // TODO + // // Imaginary part // Note: b = 0 - solver.begin_row() ; - solver.add_coefficient(u0_id, -b+d) ; - solver.add_coefficient(v0_id, -a+c) ; - solver.add_coefficient(u1_id, -d) ; - solver.add_coefficient(v1_id, -c) ; - solver.add_coefficient(v2_id, a) ; - solver.end_row() ; + // TODO + return Base::OK; } @@ -574,8 +582,10 @@ set_mesh_uv_from_system(Adaptor& mesh, // Note : 2*index --> u // 2*index + 1 --> v - NT u = solver.variable(2*index ).value() ; - NT v = solver.variable(2*index + 1).value() ; + + double u, v; + // TODO + // Fill vertex (u,v) and mark it as "parameterized" mesh.set_vertex_uv(vertexIt, Point_2(u,v)); From 6d1332c3fafdcf45ba6c8481673d0a11654014cb Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 10 Aug 2015 07:54:01 +0200 Subject: [PATCH 4/6] Replace OpenNL by Eigen Conjugate Gradient (not least squares yet). More than 2x speed up. --- Solver_interface/include/CGAL/Eigen_vector.h | 4 + .../include/CGAL/LSCM_parameterizer_3.h | 91 +++++++++---------- 2 files changed, 48 insertions(+), 47 deletions(-) diff --git a/Solver_interface/include/CGAL/Eigen_vector.h b/Solver_interface/include/CGAL/Eigen_vector.h index 6280e289be1..c8acdfe74b5 100644 --- a/Solver_interface/include/CGAL/Eigen_vector.h +++ b/Solver_interface/include/CGAL/Eigen_vector.h @@ -86,6 +86,10 @@ public: void set(std::size_t i,NT value) { this->operator[](static_cast(i))=value; } + + NT get(std::size_t i) const { + return this->operator[](static_cast(i)); + } NT* vector() {return this->data();} diff --git a/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h index 47d42c98f34..7e92a1b89b9 100644 --- a/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h @@ -24,20 +24,18 @@ #include #include - -#ifdef CGAL_EIGEN3_ENABLED -#include -#endif - #include - #include #include #include #include #include +#ifdef CGAL_EIGEN3_ENABLED +#include +#endif + #define DEBUG_TRACE /// \file LSCM_parameterizer_3.h @@ -144,12 +142,9 @@ private: typedef typename Sparse_LA::Vector Vector; typedef typename Sparse_LA::Matrix Matrix; - typedef Eigen_solver_traits::EigenType> > - LeastSquaresSolver ; - - typedef typename LeastSquaresSolver::Matrix LSS_Matrix; - typedef typename LeastSquaresSolver::Vector LSS_Vector; - typedef typename LeastSquaresSolver::NT LSS_Scalar; + // typedef typename OpenNL::LinearSolver + // LeastSquaresSolver ; + typedef EigenLeastSquaresSolver LeastSquaresSolver; // Public operations public: @@ -185,8 +180,7 @@ private: /// \pre Vertices must be indexed. /// \pre X and B must be allocated and empty. /// \pre At least 2 border vertices must be parameterized. - void initialize_system_from_mesh_border(LSS_Matrix& A, - LSS_Vector& B, + void initialize_system_from_mesh_border(LeastSquaresSolver& solver, const Adaptor& mesh) ; /// Utility for setup_triangle_relations(): @@ -198,7 +192,7 @@ private: /// Create two lines in the linear system per triangle (one for u, one for v). /// /// \pre vertices must be indexed. - Error_code setup_triangle_relations(LSS_Matrix& A, LSS_Vector& B, + Error_code setup_triangle_relations(LeastSquaresSolver& solver, const Adaptor& mesh, Facet_const_handle facet) ; @@ -301,39 +295,34 @@ parameterize(Adaptor& mesh) // Create sparse linear system "A*X = B" of size 2*nbVertices x 2*nbVertices // (in fact, we need only 2 lines per triangle x 1 column per vertex) - LSS_Matrix A (2 * nbVertices, 2 * nbVertices); - LSS_Vector B (2 * nbVertices); - + LeastSquaresSolver solver(2*nbVertices); + solver.set_least_squares(true) ; + // Initialize the "A*X = B" linear system after // (at least two) border vertices parameterization - initialize_system_from_mesh_border(A, B, mesh); + initialize_system_from_mesh_border(solver, mesh); // Fill the matrix for the other vertices + solver.begin_system() ; for (Facet_iterator facetIt = mesh.mesh_facets_begin(); facetIt != mesh.mesh_facets_end(); facetIt++) { // Create two lines in the linear system per triangle (one for u, one for v) - status = setup_triangle_relations(A, B, mesh, facetIt); + status = setup_triangle_relations(solver, mesh, facetIt); if (status != Base::OK) return status; } - - A.assemble_matrix(); - + solver.end_system() ; #ifdef DEBUG_TRACE std::cerr << " matrix filling (" << 2*mesh.count_mesh_facets() << " x " << nbVertices << "): " << timer.time() << " seconds." << std::endl; timer.reset(); #endif - + // Solve the "A*X = B" linear system in the least squares sense - LSS_Vector X (2 * nbVertices); - LSS_Scalar D; - LeastSquaresSolver solver; - if (!(solver.linear_solver (A, B, X, D))) - status = Base::ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; - + if ( ! solver.solve() ) + status = Base::ERROR_CANNOT_SOLVE_LINEAR_SYSTEM; #ifdef DEBUG_TRACE std::cerr << " solving linear system: " << timer.time() << " seconds." << std::endl; @@ -341,7 +330,7 @@ parameterize(Adaptor& mesh) #endif if (status != Base::OK) return status; - return status; + // Copy X coordinates into the (u,v) pair of each vertex set_mesh_uv_from_system(mesh, solver); #ifdef DEBUG_TRACE @@ -414,7 +403,7 @@ check_parameterize_preconditions(Adaptor& mesh) template inline void LSCM_parameterizer_3:: -initialize_system_from_mesh_border(LSS_Matrix& A, LSS_Vector& B, +initialize_system_from_mesh_border(LeastSquaresSolver& solver, const Adaptor& mesh) { for (Vertex_const_iterator it = mesh.mesh_vertices_begin(); @@ -430,13 +419,14 @@ initialize_system_from_mesh_border(LSS_Matrix& A, LSS_Vector& B, // Write (u,v) in X (meaningless if vertex is not parameterized) // Note : 2*index --> u // 2*index + 1 --> v + solver.variable(2*index ).set_value(uv.x()) ; + solver.variable(2*index + 1).set_value(uv.y()) ; // Copy (u,v) in B if vertex is parameterized - if (mesh.is_vertex_parameterized(it)) - { - B.set (2 * index, uv.x ()); - B.set (2 * index + 1, uv.y ()); - } + if (mesh.is_vertex_parameterized(it)) { + solver.variable(2*index ).lock() ; + solver.variable(2*index + 1).lock() ; + } } } @@ -491,7 +481,7 @@ template inline typename LSCM_parameterizer_3::Error_code LSCM_parameterizer_3:: -setup_triangle_relations(LSS_Matrix& A, LSS_Vector& B, +setup_triangle_relations(LeastSquaresSolver& solver, const Adaptor& mesh, Facet_const_handle facet) { @@ -554,15 +544,24 @@ setup_triangle_relations(LSS_Matrix& A, LSS_Vector& B, // // Real part // Note: b = 0 - - // TODO - + solver.begin_row() ; + solver.add_coefficient(u0_id, -a+c) ; + solver.add_coefficient(v0_id, b-d) ; + solver.add_coefficient(u1_id, -c) ; + solver.add_coefficient(v1_id, d) ; + solver.add_coefficient(u2_id, a) ; + solver.end_row() ; // // Imaginary part // Note: b = 0 + solver.begin_row() ; + solver.add_coefficient(u0_id, -b+d) ; + solver.add_coefficient(v0_id, -a+c) ; + solver.add_coefficient(u1_id, -d) ; + solver.add_coefficient(v1_id, -c) ; + solver.add_coefficient(v2_id, a) ; + solver.end_row() ; - // TODO - return Base::OK; } @@ -582,10 +581,8 @@ set_mesh_uv_from_system(Adaptor& mesh, // Note : 2*index --> u // 2*index + 1 --> v - - double u, v; - // TODO - + NT u = solver.variable(2*index ).value() ; + NT v = solver.variable(2*index + 1).value() ; // Fill vertex (u,v) and mark it as "parameterized" mesh.set_vertex_uv(vertexIt, Point_2(u,v)); From 4276d2795ecc2fe286b280866578850acac3f652 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 12 Aug 2015 11:26:53 +0200 Subject: [PATCH 5/6] Using Eigen::SimplicialLDLT. Speeds up even more (15x compared to original OpenNL impl) --- Polyhedron/demo/Polyhedron/CMakeLists.txt | 6 + .../include/CGAL/Eigen_least_squares_solver.h | 335 ++++++++++++++++++ .../include/CGAL/Eigen_solver_traits.h | 15 +- 3 files changed, 355 insertions(+), 1 deletion(-) create mode 100644 Solver_interface/include/CGAL/Eigen_least_squares_solver.h diff --git a/Polyhedron/demo/Polyhedron/CMakeLists.txt b/Polyhedron/demo/Polyhedron/CMakeLists.txt index bcf1e0e88a6..1475bae5b6f 100644 --- a/Polyhedron/demo/Polyhedron/CMakeLists.txt +++ b/Polyhedron/demo/Polyhedron/CMakeLists.txt @@ -57,6 +57,12 @@ if(Qt5_FOUND) endif(Qt5_FOUND) +find_package(OpenMP) +if (OPENMP_FOUND) + set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") +endif() + find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater) if (EIGEN3_FOUND) include( ${EIGEN3_USE_FILE} ) diff --git a/Solver_interface/include/CGAL/Eigen_least_squares_solver.h b/Solver_interface/include/CGAL/Eigen_least_squares_solver.h new file mode 100644 index 00000000000..a446c551699 --- /dev/null +++ b/Solver_interface/include/CGAL/Eigen_least_squares_solver.h @@ -0,0 +1,335 @@ +// Copyright (c) 2005-2008 Inria Loria (France). +/* + * author: Bruno Levy, INRIA, project ALICE + * website: http://www.loria.fr/~levy/software + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation, either version 3 + * of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * + * Scientific work that use this software can reference the website and + * the following publication: + * + * @INPROCEEDINGS {levy:NMDGP:05, + * AUTHOR = Bruno Levy, + * TITLE = Numerical Methods for Digital Geometry Processing, + * BOOKTITLE =Israel Korea Bi-National Conference, + * YEAR=November 2005, + * URL=http://www.loria.fr/~levy/php/article.php?pub=../publications/papers/2005/Numerics + * } + * + * Laurent Saboret 2005-2006: Changes for CGAL: + * - Added OpenNL namespace + * - DefaultLinearSolverTraits is now a model of the SparseLinearAlgebraTraits_d concept + * - Added SymmetricLinearSolverTraits + * - copied Jacobi preconditioner from Graphite 1.9 code + */ + + +#ifndef __EIGEN_LEAST_SQUARES_SOLVER__ +#define __EIGEN_LEAST_SQUARES_SOLVER__ + +#include + +#include +#include +#include + + +namespace CGAL { + + +/* + * Solves a linear system or minimizes a quadratic form. + */ +class EigenLeastSquaresSolver +{ +protected: + + enum State { + INITIAL, IN_SYSTEM, IN_ROW, CONSTRUCTED, SOLVED + } ; + +public: + // typedef Eigen_solver_traits::EigenType, + // Eigen::Lower|Eigen::Upper> > Solver; + typedef Eigen_solver_traits::EigenType> > Solver; + typedef typename Solver::Matrix Matrix ; + typedef typename Solver::Vector Vector ; + typedef typename Solver::NT CoeffType ; + + class Variable { + public: + Variable() : x_(0), index_(-1), locked_(false) { } + double value() const { return x_; } + void set_value(double x_in) { x_ = x_in ; } + void lock() { locked_ = true ; } + void unlock() { locked_ = false ; } + bool is_locked() const { return locked_ ; } + unsigned int index() const { + CGAL_assertion(index_ != -1) ; + return (unsigned int)(index_) ; + } + void set_index(unsigned int index_in) { + index_ = index_in ; + } + private: + CoeffType x_ ; + int index_ ; + bool locked_ ; + } ; + + + EigenLeastSquaresSolver(unsigned int nb_variables) { + state_ = INITIAL ; + nb_variables_ = nb_variables ; + variable_ = new Variable[nb_variables] ; + A_ = NULL ; + x_ = NULL ; + b_ = NULL ; + } + + ~EigenLeastSquaresSolver() { + delete[] variable_ ; + delete A_ ; + delete x_ ; + delete b_ ; + } + + // __________________ Parameters ________________________ + + void set_least_squares(bool x) { } + + // __________________ Access ____________________________ + + int nb_variables() const { return nb_variables_ ; } + + Variable& variable(unsigned int idx) { + CGAL_assertion(idx < nb_variables_) ; + return variable_[idx] ; + } + + const Variable& variable(unsigned int idx) const { + CGAL_assertion(idx < nb_variables_) ; + return variable_[idx] ; + } + + // _________________ Construction _______________________ + + void begin_system() { + current_row_ = 0 ; + transition(INITIAL, IN_SYSTEM) ; + // Enumerate free variables. + unsigned int index = 0 ; + for(int ii=0; ii < nb_variables() ; ii++) { + Variable& v = variable(ii) ; + if(!v.is_locked()) { + v.set_index(index) ; + index++ ; + } + } + unsigned int n = index ; + A_ = new Matrix(static_cast(n),static_cast(n)) ; + x_ = new Vector(n) ; + b_ = new Vector(n) ; + for(unsigned int i=0; iset(i, 0.); + b_->set(i, 0.); + } + variables_to_vector() ; + } + + void begin_row() { + transition(IN_SYSTEM, IN_ROW) ; + af_.clear() ; + if_.clear() ; + al_.clear() ; + xl_.clear() ; + bk_ = 0 ; + } + + void set_right_hand_side(double b) { + check_state(IN_ROW) ; + bk_ = b ; + } + + void add_coefficient(unsigned int iv, double a) { + check_state(IN_ROW) ; + Variable& v = variable(iv) ; + if(v.is_locked()) { + al_.push_back(a) ; + xl_.push_back(v.value()) ; + } else { + af_.push_back(a) ; + if_.push_back(v.index()) ; + } + } + + void normalize_row(CoeffType weight = 1) { + check_state(IN_ROW) ; + CoeffType norm = 0.0 ; + unsigned int nf = af_.size() ; + for(unsigned int i=0; i1e-40 ); + scale_row(weight / norm) ; + } + + void scale_row(CoeffType s) { + check_state(IN_ROW) ; + unsigned int nf = af_.size() ; + for(unsigned int i=0; iadd_coef (if_[i], if_[j], af_[i] * af_[j]); + } + } + CoeffType S = - bk_ ; + for(unsigned int j=0; jset(if_[i], b_->get(if_[i]) - af_[i] * S); + } + current_row_++ ; + transition(IN_ROW, IN_SYSTEM) ; + } + + void end_system() { + A_->assemble_matrix (); + transition(IN_SYSTEM, CONSTRUCTED) ; + } + + // ----------------------------- Solver ------------------------------- + + // Solves a linear system or minimizes a quadratic form. + // Return true on success. + // (modified for SparseLinearAlgebraTraits_d concept) + bool solve() + { + check_state(CONSTRUCTED) ; + + // Solve the sparse linear system "A*X = B". On success, the solution is (1/D) * X. + Solver solver_; + CoeffType D; + + bool success = solver_.linear_solver(*A_, *b_, *x_, D) ; + CGAL_assertion(D == 1.0); // WARNING: this library does not support homogeneous coordinates! + + vector_to_variables() ; + + transition(CONSTRUCTED, SOLVED) ; + + delete A_ ; A_ = NULL ; + delete b_ ; b_ = NULL ; + delete x_ ; x_ = NULL ; + + return success; + } + +protected: + + // ----------- Converting between user representation and the internal representation ----- + + void vector_to_variables() { + for(int ii=0; ii < nb_variables(); ii++) { + Variable& v = variable(ii) ; + if(!v.is_locked()) { + v.set_value( x_->get(v.index()) ) ; + } + } + } + + void variables_to_vector() { + for(int ii=0; ii < nb_variables(); ii++) { + Variable& v = variable(ii) ; + if(!v.is_locked()) { + x_->set(v.index(), v.value()); + } + } + } + + // ----------- Finite state automaton (checks that calling sequence is respected) --------- + + std::string state_to_string(State s) { + switch(s) { + case INITIAL: + return "initial" ; + case IN_SYSTEM: + return "in system" ; + case IN_ROW: + return "in row" ; + case CONSTRUCTED: + return "constructed" ; + case SOLVED: + return "solved" ; + } + // Should not go there. + CGAL_error(); + return "undefined" ; + } + + void check_state(State s) { + CGAL_USE(s); + CGAL_assertion(state_ == s) ; + } + + void transition(State from, State to) { + check_state(from) ; + state_ = to ; + } + +private: + + // --------------- user representation -------------- + unsigned int nb_variables_ ; + Variable* variable_ ; + + // --------------- construction ----------------------- + State state_ ; + unsigned int current_row_ ; + std::vector af_ ; + std::vector if_ ; + std::vector al_ ; + std::vector xl_ ; + double bk_ ; + + // --------------- internal representation --------- + Matrix* A_ ; + Vector* x_ ; + Vector* b_ ; + +} ; + + +} // namespace OpenNL + +#endif diff --git a/Solver_interface/include/CGAL/Eigen_solver_traits.h b/Solver_interface/include/CGAL/Eigen_solver_traits.h index 9e065c65558..e2d11bd09df 100644 --- a/Solver_interface/include/CGAL/Eigen_solver_traits.h +++ b/Solver_interface/include/CGAL/Eigen_solver_traits.h @@ -49,7 +49,12 @@ namespace internal { }; template - struct Get_eigen_matrix< ::Eigen::ConjugateGradient,FT>{ + struct Get_eigen_matrix< ::Eigen::ConjugateGradient,FT>{ + typedef Eigen_sparse_symmetric_matrix type; + }; + + template + struct Get_eigen_matrix< ::Eigen::SimplicialLDLT,FT>{ typedef Eigen_sparse_symmetric_matrix type; }; @@ -94,6 +99,7 @@ public: Eigen_solver_traits():m_mat(NULL), m_solver_sptr(new EigenSolverT) { + } EigenSolverT& solver() { return *m_solver_sptr; } @@ -115,6 +121,13 @@ public: X = m_solver_sptr->solve(B); + if (m_solver_sptr->info() == Eigen::NumericalIssue) + std::cerr << "Numerical Issue" << std::endl; + if (m_solver_sptr->info() == Eigen::NoConvergence) + std::cerr << "No Convergence" << std::endl; + if (m_solver_sptr->info() == Eigen::InvalidInput) + std::cerr << "Invalid Input" << std::endl; + return m_solver_sptr->info() == Eigen::Success; } From 2b676edd1e16edc7a6fe5b03d53968f39c2aa3db Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Mon, 24 Aug 2015 15:02:51 +0200 Subject: [PATCH 6/6] Clean up least square solver wrapper --- .../include/CGAL/LSCM_parameterizer_3.h | 19 +- .../internal/Eigen_least_squares_solver.h | 181 ++++++++++++++++++ 2 files changed, 189 insertions(+), 11 deletions(-) create mode 100644 Surface_mesh_parameterization/include/CGAL/internal/Eigen_least_squares_solver.h diff --git a/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h index 7e92a1b89b9..b05fde07c7f 100644 --- a/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h @@ -33,7 +33,7 @@ #include #ifdef CGAL_EIGEN3_ENABLED -#include +#include #endif #define DEBUG_TRACE @@ -142,9 +142,7 @@ private: typedef typename Sparse_LA::Vector Vector; typedef typename Sparse_LA::Matrix Matrix; - // typedef typename OpenNL::LinearSolver - // LeastSquaresSolver ; - typedef EigenLeastSquaresSolver LeastSquaresSolver; + typedef CGAL::internal::Eigen_least_squares_solver LeastSquaresSolver; // Public operations public: @@ -296,7 +294,6 @@ parameterize(Adaptor& mesh) // Create sparse linear system "A*X = B" of size 2*nbVertices x 2*nbVertices // (in fact, we need only 2 lines per triangle x 1 column per vertex) LeastSquaresSolver solver(2*nbVertices); - solver.set_least_squares(true) ; // Initialize the "A*X = B" linear system after // (at least two) border vertices parameterization @@ -419,13 +416,13 @@ initialize_system_from_mesh_border(LeastSquaresSolver& solver, // Write (u,v) in X (meaningless if vertex is not parameterized) // Note : 2*index --> u // 2*index + 1 --> v - solver.variable(2*index ).set_value(uv.x()) ; - solver.variable(2*index + 1).set_value(uv.y()) ; + solver.set_value (2*index, uv.x()); + solver.set_value (2*index + 1, uv.y()); // Copy (u,v) in B if vertex is parameterized if (mesh.is_vertex_parameterized(it)) { - solver.variable(2*index ).lock() ; - solver.variable(2*index + 1).lock() ; + solver.lock (2*index); + solver.lock (2*index + 1); } } } @@ -581,8 +578,8 @@ set_mesh_uv_from_system(Adaptor& mesh, // Note : 2*index --> u // 2*index + 1 --> v - NT u = solver.variable(2*index ).value() ; - NT v = solver.variable(2*index + 1).value() ; + NT u = solver.value (2*index); + NT v = solver.value (2*index + 1); // Fill vertex (u,v) and mark it as "parameterized" mesh.set_vertex_uv(vertexIt, Point_2(u,v)); diff --git a/Surface_mesh_parameterization/include/CGAL/internal/Eigen_least_squares_solver.h b/Surface_mesh_parameterization/include/CGAL/internal/Eigen_least_squares_solver.h new file mode 100644 index 00000000000..452a3ea1d00 --- /dev/null +++ b/Surface_mesh_parameterization/include/CGAL/internal/Eigen_least_squares_solver.h @@ -0,0 +1,181 @@ +#ifndef __EIGEN_LEAST_SQUARES_SOLVER__ +#define __EIGEN_LEAST_SQUARES_SOLVER__ + +#include +#include + +#include +#include +#include + + +namespace CGAL { + + namespace internal + { + class Eigen_least_squares_solver + { + + public: + typedef Eigen_solver_traits::EigenType> > Solver; + typedef typename Solver::Matrix Matrix ; + typedef typename Solver::Vector Vector ; + typedef typename Solver::NT Scalar ; + typedef typename CGAL::cpp11::tuple Indexed_scalar; + + + private: + + std::vector _indexed_scalars; + + std::vector _af; + std::vector _if; + std::vector _al; + std::vector _xl; + + Matrix* _A ; + Vector* _x ; + Vector* _b ; + + public: + + Eigen_least_squares_solver (unsigned int n_bvariables) + { + _indexed_scalars.resize (n_bvariables); + for (std::size_t i = 0; i < _indexed_scalars.size (); ++ i) + { + _indexed_scalars[i].get<0>() = 0.; + _indexed_scalars[i].get<1>() = -1; + _indexed_scalars[i].get<2>() = false; + } + _A = NULL ; + _x = NULL ; + _b = NULL ; + } + + ~Eigen_least_squares_solver() + { + if (_A != NULL) delete _A ; + if (_x != NULL) delete _x ; + if (_b != NULL) delete _b ; + } + + const Scalar& value (std::size_t index) const + { + return _indexed_scalars[index].get<0>(); + } + + void set_value (std::size_t index, Scalar value) + { + _indexed_scalars[index].get<0>() = value; + } + + void lock (std::size_t index) + { + _indexed_scalars[index].get<2>() = true; + } + + void store_scalars_in_x () + { + for (std::size_t i = 0; i < _indexed_scalars.size (); ++ i) + if (!(_indexed_scalars[i].get<2>())) + _x->set (_indexed_scalars[i].get<1>(), _indexed_scalars[i].get<0>()); + } + + void get_result () + { + for (std::size_t i = 0; i < _indexed_scalars.size (); ++ i) + if (!(_indexed_scalars[i].get<2>())) + _indexed_scalars[i].get<0>() = _x->get(_indexed_scalars[i].get<1>()); + } + + void begin_system () + { + std::size_t index = 0; + for (std::size_t i = 0; i < _indexed_scalars.size (); ++ i) + if(!(_indexed_scalars[i].get<2>())) + _indexed_scalars[i].get<1>() = index ++; + + _A = new Matrix (index, index); + _x = new Vector (index); + _b = new Vector (index); + + for (std::size_t i = 0; i < index; ++ i) + { + _x->set(i, 0.); + _b->set(i, 0.); + } + + store_scalars_in_x() ; + } + + void begin_row () + { + _af.clear() ; + _if.clear() ; + _al.clear() ; + _xl.clear() ; + } + + void add_coefficient (unsigned int iv, double a) + { + if(_indexed_scalars[iv].get<2>()) + { + _al.push_back(a) ; + _xl.push_back(_indexed_scalars[iv].get<0>()); + } + else + { + _af.push_back(a) ; + _if.push_back(_indexed_scalars[iv].get<1>()); + } + } + + void end_row () + { + std::size_t nf = _af.size() ; + std::size_t nl = _al.size() ; + + for (std::size_t i = 0; i < nf; ++ i) + for (std::size_t j = 0; j < nf; ++ j) + _A->add_coef (_if[i], _if[j], _af[i] * _af[j]); + + Scalar S = 0.; + + for (std::size_t j = 0; j < nl; ++ j) + S += _al[j] * _xl[j] ; + + for (std::size_t i = 0; i < nf; ++ i) + _b->set(_if[i], _b->get(_if[i]) - _af[i] * S); + } + + void end_system () + { + _A->assemble_matrix (); + } + + bool solve () + { + Scalar D; + Solver solver; + bool success = solver.linear_solver(*_A, *_b, *_x, D) ; + + get_result() ; + + delete _A ; _A = NULL ; + delete _b ; _b = NULL ; + delete _x ; _x = NULL ; + + return success; + } + + + + } ; + + + } // namespace internal + +} // namespace CGAL + +#endif